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ABSTRACT 

Empirical "backward" galaxy evolution models for infrared bright galaxies are constrained using 
multi-band infrared surveys. We developed a new Monte-Carlo algorithm for this task, implementing 
luminosity dependent distribution functions for the galaxies' infrared spectral energy distributions 
(SEDs) and for the AGN contribution, allowing for evolution of these quantities. The adopted SEDs 
take into account the contributions of both starbursts and AGN to the infrared emission, for the first 
time in a coherent treatment rather than invoking separate AGN and star-forming populations. 

In the first part of the paper we consider the quantification of the AGN contribution for local 
universe galaxies, as a function of total infrared luminosity. It is made using a large sample of LIRGs 
and ULIRGs for which mid-infrared spectra are available in the Spitzer archive. We find the ratio 
of AGN 6 /im luminosity and total infrared luminosity to rise with i^iR over the infrared luminosity 
range 10^^ to 10^^ Lq and estimate its spread. Judging from the modest number of distant sources 
with Spitzer spectroscopy, the relation changes at high z. 

In the second part we present the model. Our best-fit model adopts very strong luminosity evolution, 
L = Lo{l + z)^ "*, up to z = 2.3, and density evolution, p = pq{1 + z)^, up to z = 1, for the population 
of infrared galaxies. At higher z, the evolution rates drop as (l-fz)"^ and respectively. To 

reproduce mid-infrared to submillimeter number counts and redshift distributions, it is necessary to 
introduce both an evolution in the AGN contribution and an evolution in the luminosity-temperature 
relation. At a given total infrared luminosity, high redshift infrared galaxies have typically smaller 
AGN contributions to the rest frame mid-infrared, and colder far-infrared dust temperatures than 
locally. We also suggest an extension of the local infrared galaxy population towards lower dust 
temperatures. 

Our models are in plausible agreement with current photometry-based estimates of the typical AGN 
contribution as a function of mid-infrared flux, and well placed to be compared to upcoming Spitzer 
spectroscopic results. As an example of future applications, we use our best-fitting model to make 
predictions for surveys with Herschel. 

Subject headings: cosmology: source counts, cosmology: redshift distribution, galaxies: active, galax- 
ies: starburst, galaxies: evolution, infrared: galaxies 



L INTRODUCTION 

The disc overy of the cosmic i nfrared background 
(GIB) (see iHauser fc Dwkl [200ll for a review) , to- 
gether with recent deep cosmological surveys in the 
infrared (IR) and submillimeter bands, have offered 
new perspectives on our understanding of galaxy for- 
mation and evolution. The surprisingly large amount 
of energy contained in the GIB showed that it is cru- 
cial to probe its contributing galaxies to understand 
when and how the bulk of stars formed in the Uni- 
Thanks t o the deep co smological surveys car- 
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Spitzer dW erncr ct al. " 2004 e.g iPapovich et all 12004 
iFraver et al.. .2006a: Dole et al.l |2006[ ) it is now pos- 
sible, to various degrees, to resolve the GIB into 
discrete sources. The source counts are high 
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when compared to no-evolutio n or m oderate-evolution 
models for infrared ga laxies (jCuid crdoni ct al. 199^; 
iFranceschini et al]ll998f ). The striking results of these 
surveys concerning the evolution of the infrared and sub- 
millimeter galaxy population, and the constraints from 
the measurements of the GIB, require the development 
of new models to explain the high rate of evolution of 
infrared galaxies. 

The problem of describing the number distribution of 
galaxies in the Universe is usually tackled via one of two 
methods. In the first method, known as forward evolu- 
tion, the calculation assumes some initial conditions and 
physical processes for chemical and photometric evolu- 
tion of the stellar populations that heat the dust. The so- 
called semi-analytical approach combines these assump- 
tions about the chemical/photometric evolution of galax- 
ies with models for the dissipative and non-dissipative 
processes affecting galaxy formation within dark matter 
halos (Baugh et al..i200& ) and has pro vided a reasonable 
fit to the source counts in the infrared (jGuiderdoni et al.l 
ll998l : [Lacev et al.ll2008[ ). Arguably, the forward evolution 
approach has the advantage of being based on a more 
fundamental set of assumptions. However the obvious 



lisadvantage is the large number of free/unknown pa- 
'ameters assumed in these models. Compared to other 
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wavelengths, there is one further comphcation for for- 
ward modeUing of infrared surveys. Even if global prop- 
erties of the evolving galaxies like star formation rates 
or AGN activity are correctly modelled, additional as- 
sumptions about the dust content and structure of the 
galaxy have to be invoked to convert them into a predic- 
tion of the luminosity and spectral energy distribution of 
the re-emitted infrared radiation. 

The alternative method, often called the backward evo- 
lution approach, takes the observed, present day (z = 0) 
luminosity function (LF) and evolves it in luminosity 
and/or density out to higher redshi fts assuming some 
param e trization of the evolution fe.g. lFranceschini at al 
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1988t JP earson fc Rowan-Robinson 

2001; Franceschini et all 1200 ll : iRowan-RobinsonI l200ll 
Lagache et al. 200^. This method has the advantage of 
being both direct and relatively simple to implement. 
The disadvantage in the past was that the information 
on which the measurements of evolution were often 
made with IRAS or ISO data which extended out only 
to relatively low redshifts. Because of their simplicity, 
backward evolution models have traditionally played a 
strong role in the planning of new infrared surveys. A 
backward evolution model fitting the main constraints 
provided by previous missions can be easily modified to 
predict the results of new missions and help in the first 
steps of interpreting their results. Thanks to the latest 
deep and comprehensive observations cited above, the 
tools to constrain the backward evolution methodology 
to significantly higher redshifts are now available. 

Very recently, several empirical approaches have been 
proposed to model the high rate of evolution of infrared 
galaxies, in particular to reproduce sourc e counts of the 



mid-IR surve ys made with Sp itzer Ce .g. iLagache et al 
20031 120041: iGruppioni et al.l 120051 : iRowan-Robinson 



20091) . 

The "classical" backward evolution model starts from 
different populations of galaxies, typically cirrus, star- 
burst, active galactic nuclei (A GN) and ultraluminou s 
infrared galaxies (ULIRGs) (e.g. lRowan- Robinsonl[200l . 
or a subset of those. Each population is assigned a spec- 
tral energy distribution (SED) and local luminosity func- 
tion, and evolved independently, assuming it has a proper 
cosmic evolution rate. These hypotheses are not always 
satisfied. First of all, the SEDs adopted for starburst 
galaxies are usually either represented by a single SED or 
come from a set of templates where a unique relation be- 
tween temperature and lumi nosity is assumed. Thi s as- 
pect was already dis cussed bv lChapman et all (|2003[ ) and 
IChapin et al.l ()2009l ): analyzing a sample of local infrared 
galaxies, they found that the luminosity-temperature 
{L — T) relation presents a significant spread. The lat- 
ter, in particular, showed evidence for some evolution 
of the relation with redshift, comparing the local trend 
with high-redshift dat a from SHARC-II 350 Mm obser- 
vations of SMGs (e.g. iKovacs et al.ll2006t ICoppin et al.l 
[200 81 . Second, and perhaps more important, different 
populations of galaxies are not as distinct as assumed in 
the models: very often AGN and starbursts co-exist in 
the same object, can be predominant at different times 
depending of the evolutionary stage of the object itself 
and can influence each other (e.g. AGN feedback on star 
formation) . 

Backward evolution models based on single SEDs or 



the simple SED family approach described above have 
been quite successful in fitting the number counts from 
the IRAS and ISO missions, the first submillimeter 
counts, the global GIB level, and for making predictions 
for Spitzer and Herschel. The number and quality of ob- 
served constraints is increasing, however, and already the 
first Spitzer results have led to on-t he-spot modifications 
(c.f. the modified SEDs adopted bv ILagache et al.|[2004D 
that may either genuinely represent improved knowledge, 
or reflect limitations of the simple assumptions made pre- 
viously. Furthermore, questions gaining increased impor- 
tance, like the co-existence of AGN and star-formation 
in infrared galaxies, cannot be addressed by this gener- 
ation of models, not even in the simple sense of fitting 
existing data and extrapolating to new observations. The 
goal of this work is hence to take the next step and de- 
velop a backward evolution model that considers realistic 
spreads in far-infrared SEDs, in AGN contributions at 
different luminosities, and their possible evolution with 
redshift. To that end, we first have to describe the lo- 
cal situation that is the starting point of the backward 
evolution scenario. 

Luminous (LIRGs: 10" < Ltir = 
10^2 Lq) and ultr aluminous (ULIRGs : Ltir > 10^^ Lq) 
infrared galaxies ([Sanders et al.l ll988^ have been studied 
extensively in the local Universe, for insta nce with the 
Infra r ed Astronomical Sate llite (IRAS; e.g. ISoifer et al.l 
Il987t [Saunders et all 119901), the Infrared Space Obse r- 



vator y (ISO; e.g . Lutz et all [1998[ : [Genzel fc Gesarskvl 
[20001 : iTran et al.l[2001f ), a nd more recently, w ith the In- 
frared Spectrograph (IRS; i Houck et al .l |200_4D on board 



of Spitzer (e.g. 


Wccdman et al.l 


2005: Brandl et al. 
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2007>:.Veilleux et al. 


20091). These galaxies exhibit a large range of properties 



in the mid-IR, some showing strong PAH emission fea- 
tures characteristic of p owerful (up to 1000 Mf7)yr~^) 
star f ormation rates (e.g. lBrandl et al.ll2006l : [Smith et al.1 
l2007t l. and exhibiting a large rang e in 9.7 Mm silicate ab- 
sorption or e missi on strengths (e.g.lWeedman et al.ll2005[ : 
iDesai et al.l [20071 : [Imanishi et al 2007^. Spitzer IRS is 
now enabling the study of mid-IR spectr a of infrared 
galax i es to much higher redshifts (z > 2; e.g.lHouck et al.1 
20051: lYan et aLl [20071: [Menendez-Delmestre et al.l [20071: 



Valiante et all l2007t fPope et all l2008f K Although lo- 
cally rare, infrared galaxies are an important popula- 
tion at high redshifts and account for an increasing 
fraction of the star-f ormation activity in the universe 
([Le Floc'h et al.l[2005[ ) . By studying their infrared prop- 
erties, we are just starting to estimate the extent to which 
AGN and star-formation contribute to their infrared lu- 
minosities, and therefore determine a correct census of 
starbursts and AGN at epochs in the Universe when their 
luminosity density was at its maximum. 

The local luminosity function of IR-bright galaxies 
can be described considering two different populations: 
normal galaxies, dominating the "low-luminosity" part 
of LF, and starburst galaxies, dominant in the "high- 
luminosity" part. The AGN contribution appears domi- 
nant only at very high luminosity (Ltir ^ 2 x IO^^Lq). 
Nevertheless, a small fraction of the total infrared lu- 
minosity can be due to the presence of an active galac- 
tic nucleus even in star forma tion dominated cases, and 
vice versa ([Genzel et al.|[l998[) . There are several stud- 
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ies of the AGN content for luminous galaxies at ULIRG 
and HYLIRG (hyperluminous infrare d galaxies, Ltir. > 
10^'^ L p,) levels , both using X-ray (iFranceschini et aT 



2003; 



19981: 



Teng et al . 2005) and mid-IR (e. g. iGenzel et al 



Lutz et all 11998: Tra n et all l2001h emission, and 



it is generally believed that AGN are typically less im- 
portant at lower infrared luminosities. Still missing 
is a comprehensive study including lower luminosities 
(< 10^^ Lq), with the aim of quantifying the distribution 
of AGN luminosity respect to the infrared luminosity of 
the host and the fraction of infrared luminosity due to 
accretion. Such a study is made in the following section 
(see §[2]). 

This work is organized as follows: in the first part 
(§[2]), starting from an IRAS-based sample of LIRGs and 
ULIRGs observed with IRS, we derive local distribution 
functions for the AGN contribution at different infrared 
luminosities. Using a small sample of distant galaxies, we 
also try to explore this relation at higher z. In the second 
part (§[3|), we present a new backward evolution model 
whose new element is to treat infrared galaxies as a single 
population with a single local luminosity function, but 
with realistic spreads in both SED properties and AGN 
content. An SED including the contribution due to star- 
bursts and AGN is associated with each source, following 
the relations derived in the first part for local and distant 
objects. This model, besides reproducing existing source 
counts, redshift distributions and GIB intensity, is also 
able for the first time to quantify the contribution due to 
starbursts and/or AGN to the total infrared luminosity 
and how this contribution evolves for different luminosity 
classes. 

2. THE AGN CONTRIBUTION IN LOCAL ULIRGS AND 
LIRGS 

In order to be useful for the backward modelling 
project outlined above, a local calibration of the dis- 
tribution of AGN contributions at given infrared lumi- 
nosity has to span the highest luminosities down to at 
least the low end of the LIRG regime (lO^^i©). Given 
the evidence for the major contribution of objects above 
this luminosity threshold to the GIB and to the cos- 
mic star formation rate at z > 1 (e.g. lElbaz et al.l 
[200llPerez-Gonzalez et all[2nn5h . covering this luminos- 
ity range is essential. Historically, many of the references 
cited above have focusscd on the ULIRG (> lO^^i©) 
regime and even its upper end where the AGN con- 
tribution is highest, while studies at lower luminosities 
mostly focussed on individual interesting objects. What 
is needed is a quantification of AGN content for an unbi- 
ased far-IR selected sample reaching down to 10^^ Lq. Of 
the two principle routes towards this observational goal. 
X-ray observations and mid-IR spectroscopy, we make 
use of the second for two reasons. First, there is not yet 
a full unbiased and deep X-ray dataset for such a sam- 
ple available. There are XMM and Ghandra data of the 
depth required for good X-ray spectral analysis for many 
local LIRGs/ULIRGs, but still with a tendency to tar- 
get known AGN. Second, since the implicit goal of our 
modelling effort is to characterize the AGN effect on the 
infrared SED, mid-IR data can provide useful constraints 
quite directly even when adopting simple analysis meth- 
ods, while evidence from other wavelengths would have 
to go through a conversion to the mid-IR range, or even 



in two steps via the AGN bolometric power, a process 
that will be subject to the considerable scatter of AGN 
SEDs. 

Our basic approach is to determine the distribution of 
AGN contributions for local LIRGs and ULIRGs from 
the I RAS Revised Bright Galaxy Sample (Sanders et aP 
I2003D for which archival Spitzer spectroscopic coverage 
is essentially complete. We adopt a simple quantification 
of the AGN contribution from the rest frame 6 fim con- 
tinuum which, while not making use of the full details of 
the best S/N spectra, is applicable to the entire sample 
and allows for straightforward use in the later modelling. 

In mid-IR spectra of AGN, the 6 /im flux after subtrac- 
tion of the starbursts contribution (/6AGn), is related to 
the nuclear activity, as indicated by the intrinsic hard 
X-ray fiux ijLutz et al.l[2004 ). Hard X-rays, unless ex- 
tremely obscured in fully Gompton-thick objects, can 
provide a direct view to the central activity, while the 
infrared continuum is due to AGN emission reprocessed 
by dust, either in the torus or on somewhat larger scales. 
The mid-IR continuum has the advantage of showing no 
significant differences between type 1 and type 2 AGN 
and of being a good tracer of nuclear activity even in 
those cases where hard X-rays are strongly absorbed. 
Low resolution mid-IR spectra of galaxies can be de- 
composed into three components (Laurent et al. 20Q3): a 
component dominated by the aromatic "PAH" features 
arising from photodissociation regions or from the dif- 
fuse interstellar medium of the host, a continuum rising 
steeply toward wavelengths beyond 10 /im due to HII re- 
gions and a typically fiatter thermal AGN dust contin- 
uum present in active galaxies. 

In the range covered by the low resolution modules of 
IRS (SLl and SL 2, from ~5 to ^lA^m as described in 
iHouck et al.ll200^ . the AGN emission is most easily iso- 
lated shortwards of t he complex of aromatic emission fea- 
tures (jLaurent et al . 2000) . We determine the continuum 
at 6 /xm rest wavelength and eliminate non-AGN emission 
by subtracting a star formation template scaled with the 
strength of the aromatic "PAH" features arising from the 
host or from circumnuclear star formation. This method 
does not require spatially resolving the AGN from the 
host and, thanks to the high sensitivity of IRS, allows 
the detection of a weak AGN in the presence of strong 
star formation in many nearby galaxies. 

For the purposes of this work, the use of mid-IR spec- 
tra, in particular of /eAGN, is the best way to quantify the 
AGN fraction for several reasons: (1) in order to derive 
a local luminosity function, we need a large sample: the 
number of objects for which mid-IR spectra are available 
is much greater than the number of sources with X-ray 
data of sufficient quality; (2) we are interested in the 
contribution of AGN to the infrared emission: -^^agn is 
more strongly correlated with the AGN's infrared emis- 
sion than, for example, L2-iokcv; (3) /oagn is the only 
consistent way to constrain the AGN contribution in our 
comparison sample at high redshifts, due to the limi- 
tations (in wavelength coverage and S/N) of the high-z 
spectra; (4) as explained above, /eAGN is also sensitive to 
Gompton-thick AGN, not detected in the X-rays: these 
kinds of galaxies become in creasingly importa nt at high- 
z, as shown both in mid-IR (jDaddi et al.ll2007h and X-ray 
(|Gim et al.ii2007 ) surveys. 
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2.1. Sample selection and data reduction 

The starting point for this study is the IRAS Re- 
vised Bright Galaxy Sample (RBGS), a complete flux- 
limited survey of all extragalactic objec ts with 60 /im 
flux /so > 5.24 Jy (jSanders et all 120031) . The galax- 



ies with Ltir > 10"'^^ Lq were selected from this sam- 
ple and the mid-IR IRS data were obtained from the 
Spitzer archive for all of the them, with the exception of 
6 sources that were not yet observed or for which there 
were no observation programs at all. Most of our ob- 
servations have been obtained within the program PID 
30323 (PI L.Armus). A comprehensive study of LIRGs 
and ULIRGs is being carried on by the Great Observa- 
tory All-sky LIRG Survey (GOALS'*) and a more com- 
plete approach to derive AGN contributions from mid-IR 
spectra will be presented in forthcoming studies by the 
GOALS team. Additional data have been retrieved from 
other Spitzer programs'^. 

Because of the limit in flux, this sample misses the 
most luminous infrared galaxies that are rare in the lo- 
cal universe. In fact, it does not include any galaxy with 
L > 1O^^-^L0, excluding all HYLIRGs, for example, from 
the analysis. In order to enlarge our range in luminos- 
ity, in particular including representatives of the most 
luminous tail of the local galaxy population, we added 
Spitzer spe c tra fo r 14 galaxies from the sample of 16 by 
iTran et all (|200lD . These latter galaxies were selected 
from surveys at different limits in /go and with a prefer- 
ence for high- luminosity targets. They do not introduce 
any bias in mid-IR spectral properties or AGN content. 
In particular, no AGN-related IRAS color criteria, like 
the /25//6O ratio, were applied in the selection. 

The total sample includes 211 sources spanning a lu- 
minosity range 10*°-^ < £tir < 10*^-^ Lq (see Fig.[T]and 
Tab. ??). 

We reduced the data as follows. We subtracted, for 
each cycle, the two nod positions of the basic calibrated 
data (BCD) frames. In the difference just calculated, 
we replaced deviant pixels by values representative of 
their spectral neighbourhoods. We subtracted residual 
wavelength dependent backgrounds, measured in source- 
free regions of the two dimensional difference spectra. 
In averaging all the cycles of the 2-dimensional sub- 
tracted frames, we excluded values more than three times 
the local noise away from the mean. The calibrated 1- 
dimensional spectra for the positive and the negative 
beams were extracted and the two 1-dimensional spec- 
tra averaged in order to obtain the final spectrum, using 
the SMART analysis package (Higdon et al. 2004). 

We derive /gagn by a decomposition over a range from 
5.5 /Ltm to 6.85 fim rest wavelength. We fit the spectrum 
by the superposition of a star formation component dom- 
inated by the 6.2 fim PAH feature (M82 spectrum by 
ISturm et al.ll206of) and a simple linear approximation of 
the AGN continuum. Fig. [2] shows the decomposition of 
the spectrum for sources with different nature: NGC1614 
is a pure starburst; NGC1365 is classified as a Seyfert 
1.8 from its optical spectrum, h as a Compton-thin A GN 
observed in the X-ray emission (jRisaliti et al.l[2000l ) but 

* http://goals.ipac.caltech.odu/ 

5 PIDs 14, 105, 666 (PI J.R.Houck), PID 1096 (PI H. Spoon), 
PID 3237 (PI E.Sturm), PID 3247 (PI C. Struck), PID 20549 (PI 
R.Joseph), PID 20589 (PI C.Leitherer) 



shows also PAH features in the mid-IR spectrum; MCG- 
03-34-064 has a typical AGN-dominated mid-IR spec- 
trum; the galaxy pair NGC3690/IC694 (Arp299) is clas- 
sified as Compton-thic k AGN from its X-ray spectrum 
([Delia Ceca et al1l2002l ) and its mid-IR spectrum reveals 
that the nuclear activity is quite high but coincident with 
strong star formation. 

The results of the fitting procedure are /eAGN, repre- 
sented in Fig. [2]by the thick vertical line, and /6.2pcakSB, 
an average flux density of the PAH component over the 
rest wavelength range 6.1 to 6.35 /im (see Tab. ??). A 
more comprehensive and in-depth study of the PAH fea- 
tures of this sample will be presented in future papers by 
the GOALS team. The spectrum of M82 is technically 
well suited for this decomposition but represents just 
a single object that may not be representative for star 
forming objects in gene ral. We dea l with this using the 
results and approach of iLutz et al.l (|2004l ) (their section 
2.1). Specifically, we first use the M82 SWS spectrum 
for decomposition and then correct for systematic differ- 
ences between M82 and a sample of 11 starbursts. To 
take into account the dispersion among star-forming ob- 
jects in r elative imp ortance of al l the com ponents in the 
model of I Laurent et al. (200l), iLutz et a l. (2004) have 
decomposed this small sample of 11 star- forming objects 
using the same approach and same M82 template. They 
find a mean residual of 0.096 x /6.2pcakSB with a disper- 
sion of 0.085. This corresponds to a a 6 fim continuum in 
M82 that is fainter than in other starbursts. T his finding 
agrees with the results of iBrandl et al.l (|2006D who find 
more 6 /J,m continuum in their Spitzer starburst template 
composed of 13 source s than in the M82 S WS spectrum. 
Decomposition of the iBrandl et al.l ()2006l ) template us- 
ing our approach and correction for M82 systematics in 
the way described results in an insignificant residual con- 
tinuum of less than 0.5% of the 6.2fim PAH peak, show- 
ing that also quantitatively the agreement between their 
template and our approach is good. Tab. ?? implements 
our approach in the following way: the values /gagn 
used for all further analysis are obtained after subtract- 
ing 0.096 X /6.2pcakSB from the direct fit result /eACN- 
They are thus corrected for the slight systematic differ- 
ence between M82 and other starbursts. Error estimates 
for /eAGN are the quadratic sum of two components. The 
first is a measurement error based on individual pixel 
noise derived from the dispersion in the difference of ob- 
servation and fit. The second is 0.085 x /6.2pcakSB, thus 
considering the dispersion in the properties of the com- 
parison star forming galaxies. 

In summary, AGN 6 /zm continua obtained by our de- 
compositi on, are thus c onsistent with those obtained us- 
ing Brand Tet al.l (j2006[ ) starburst template, and in addi- 
tion include an estimate of the uncertainty arising from 
the spread of starburst properties. 

2.2. Distribution of the AGN contribution 

Next we study the AGN contribution to Ltir in 
ULIRGs and LIRGs as a function of luminosity. Our 
results will be used in the modeling of § [3] to build com- 
posite SEDs taking into account the two main energy 
sources (starbursts and AGN) of infrared galaxies. 

Fig. [3] shows the distributions of the ratio between the 
6 fim AGN continuum luminosity and the total infrared 
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luminosity (vLqagn / L^ib,) for five different luminosity 
bins, spanning 10^°'^ to 10^^'^ Lq. These distributions do 
not show the intrinsic AGN luminosity, but the fraction 
of the IR emission due to the contribution of an AGN, 
because of the normalization used. 

As shown from the distributions of Fig. [3] (solid and 
cross-hatched histograms) and also in the diagram of 
Fig. [7] (open squares), there is no evident trend with lu- 
minosity in the values of the means of the detections. In 
the bins at lower Xtir there are good statistics, but the 
percentage of upper limits to LgAGN is rather high (up- 
per limits are represented by cross-hatched histograms). 
In the bins at higher Ltir, where the /bagn detection 
rates increase in a significant way (solid histograms rep- 
resent detections), the numbers of objects are more mod- 
est. In both cases, the averages obtained from detections 
only (open squares in Fig. [7]) have to be considered up- 
per limits to the real mean values of the distributions of 

The clearest trend in our measurements is the detec- 
tion rate as a function of Ltir- The decreasing frac- 
tion of upper limits clearly demonstrates that the distri- 
bution is changing across the luminosity bins and that 
there must be a trend also in the intrinsic distribution 
of vLexGN / Ltiy{ . With the assumption of a normal dis- 
tribution, it is possible to estimate the main parameters, 
mean and standard deviation, of the intrinsic distribu- 
tion of vLqagn / Ltib by Monte-Carlo simulations, tak- 
ing into account the flux limit defining the RBGS sample 
as well as the Spitzer-IRS detection limit. 

2.3. Simulating the distribution of the AGN 
contribution 

In order to reproduce our measurements, in terms of 
detection rate, mean and standard deviation of the detec- 
tions (hereafter d.r., rridot, Cdot), using Monte-Carlo sim- 
ulations, we need to make assumptions about the intrin- 
sic distribution of uLqagn/Ltir we are looking for and 
the quantities that affect our upper limits. uLqagn/Ltir 
is tightly correlated with two different quantities. Lqagn 
depends in a trivial way on /eAGN, i-e. on the sensitivity 
of the instrument. Ltir is related to the 60 fim flu x by re- 
lations well known and tested on IRAS data (jHelou et al.l 
[19851: iDale eFalllMilh . 

Our RBGS-based sample is selected at 60 /xm. There- 
fore, the flrst step of the simulation consists of generat- 
ing a population of sources having a 60 fim flux density 
spanning a ran ge from 5.24 to 200 Jy and a distribution 
NiU) oc (jSanders et aI]|2003D . expected for a com- 

plete sample of objects in a non-evolving Euclidean uni- 
verse that is a reasonable approximation for the relatively 
small redshift range covered by our sample. We derive 
^TiR for each source generated starting from its 60 /xm 
flux density, adopting Ltir ~ 2 x J^-Leo^m as appropri- 
ate for our sample. Assuming that the distribution of 
i/LgAGN /^TiR is gaussian with a certain mean and stan- 
dard deviation (m, cr), we calculate the expected value of 
/eAGN- At this point we apply a detection limit to /gagn 
and calculate d.r., mdot and (Tdct- The values of m and a 
that best reproduce the measured d.r., mdct and Cdot for 
each luminosity bin are assumed to be the parameters of 
the i^igAGNZ-^TiR intrinsic distribution. 

Our IRS data come from the public Spitzer archive, 
consisting of several observing campaigns to varying 



depths. To assign a detection limit to /eAGN, we con- 
sider the values of 3a for the detections and the values of 
the upper limits for the non-detections in each luminosity 
bin. We assume the median values of the distributions of 
these quantities as the detection limits of our measure- 
ments. These values span ^18 — 1 mJy, depending on 
the luminosity bin (see Tab.[T|). 

Fig. |4] shows the results of the simulation for a popu- 
lation of 2 X 10'' objects. This number was chosen such 
that the sampled distributions converge. Each diagram 
shows the best gaussians that, given the detection limit, 
reproduce the observed values of d.r., rridct and CTdct- A 
summary of the distribution parameters of measurements 
and simulations is presented in Tab.[TJ The values found 
for m and a are also plotted in Fig. [7] (green stars). Here, 
a correlation between the contribution to the infrared 
luminosity due to the AGN and the total infrared lu- 
minosity is clearly visible, with a dispersion decreasing 
with the luminosity. The best fit, uLqagn / Ltir oc iTiRi 
gives a = 1.4 ± 0.6. 

Our results assume the i^Lsagn/^tir distribution to 
be gaussian for each luminosity bin. This is a reason- 
able flrst hypothesis and it reproduces the observations 
in terms of detection rate, mean and dispersion of the 
detections with a small numbers of parameters. More 
attention is needed when applying it to our Monte Carlo 
simulations of the infrared sky. The low AGN activ- 
ity end of the distribution, corresponding mostly to IRS 
non-detections, is certainly poorly constrained, but such 
differences between already weak and totally insignifi- 
cant AGN will not have a major effect on the predictions 
of our model. Of course, using these distributions to 
study physics of low luminosity AGN in these systems 
would be misleading. Differences at the high AGN activ- 
ity end can be more important by introducing or miss- 
ing very luminous mid-IR sources dominated by AGN. 
Comparing for example the log(i^L6AGN/iTiR) distribu- 
tion for the third luminosity bin in Fig. [3] and Fig. [?] 
(11.7 < log(LTiR) < 12.1), we find the modelled distri- 
bution extending up to log(i>LQAGN / L^ib.) ^ which is 
not reached by the observations. This could reflect lim- 
its in the statistics of our IRS sample as well as a true 
overprediction of our simple gaussian hypothesis at the 
high log(i'i6AGN/-^TiR) end. We will return to this issue 
in i) 14.21 in the context of our backward evolution model. 
In the following we assume that the increase of the frac- 
tion of AGN occurs only up to the last luminosity bin 
measured (Ltir ~ 10^^ i©) and then the relation be- 
tween z/igAGN / Ltir and -Ltir becomes flat (Eq. [1] see 
also Fig. [71 green dotted line) . 

2.4. Evolution with redshift 

It is interesting to check if and how the relation just 
found between vLqagn / ^tir and Ltir evolves with red- 
shift. Previous studies already showed that, even if in 
the local universe ULI RGs tend to be AGN-dominated 
at L tir > lO^'^ ;^ L^r, jLutz et al.l 119981 : IVeilleux et all 
119991 iTran et al.ll2001[ ). this trend does not extend to 
higher redshifts. For example, submillimeter galaxies 
that, with their typical infrared luminosity of 10^^ Lq, 
are among the most luminous objects known, are mainly 
starburst dominated. Mid-IR spectroscopy in fact shows 
that the contribution to the total luminosity due to 
AGN continuum is small in most of the sources ob- 
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served with Spitzer IRS at high redsh i fts (iLutz et al 



2005at iMenendez-Delmestre et alJ l2007t IVahante etal 
20071: iPope et all 120081 1 These results are consis 



tent with X-ray obser vations of the same population 
([Alexander et al.l l2005V 

A systematic study at high redshifts is difhcult with ex- 
isting data. Ideally, to be consistent with the locally de- 
rived IRAS luminosity function, a high-z sample should 
be selected at 60 x (1 -|- z) //m, at the peak of the dust 
emission, but this cannot be done with our pre-Herschel 
observations. 

Our high reds hift sample includes all th e source s 
observed by iMenendez-Delmestre et all (|2007l ). 
Vahante et all (l2007l) . iPope et all (|2008[ ) and 
Brand et al.l ( 20081 ) for which the observations cover 
5.5 — 6.85 /zm in the rest frame (see Tab. [T]). Objects 
from the first three works are submillimeter galaxies, 
some of them radio pre-selected. These galaxies are 
mainly starburst-dominated, but can present continuum 
emission in their mid-IR spectra, even if AGN is not 
the dominant power source. Galaxies from ,Brand et al.l 
([^OS) are optically faint objects selected at 70 ^m: 
their mid-IR spectra show that they can be either PAH 
or absorption dominated. All the IRS data have been 
retrieved from the Spitzer archive and reduced using the 
same techniques described in § I2.H in order to derive 
-^6AGN in a consistent way for all the objects. The 
results of the fit process are listed in Tab. [T] The values 
of Ltir have been taken from the works cited above. 

Fig. [5] shows the distributions of the vLqagn / Ltib.- 
Unfortunately, the statistics are rather modest for the 
bins at lowest and highest luminosities. The redshift 
range spanned by this sample is very broad (0.37 < 
z < 3.35). In particular, the bin at lowest luminosity 
is mainly populated by objects lying at z < 1, while the 
most luminous bins include the most distant sources. We 
treat the whole sample as a unique population, because 
the inclusion of the z ^ 1 objects does not significantly 
change the vLqagn / Ltiu distribution in any of the lu- 
minosity bins. 

In order to characterize the intrinsic distribution of 
i^LsagnZ-^tir, we run Monte-Carlo simulations with the 
same techniques used for the local galaxies explained in 
§ 12.31 We generate a population of sources having 60 ^m 
flux density spanning 40 to 1100 mJy. We calculate the 
expected value of the 60 /im rest-frame flux density from 
the 850 /im and 70 fim fluxes, respectively for the sub- 
millimeter and the 70 /zm selected sources, assuming a 
grey-body spectrum with emissivity index /3 = 1.5 and 
Td — 32K. The values obtained give us the range of the 
60 ^m population to simulate. Because of the low counts 
in the bin at highest luminosity, we are able to repro- 
duce only the distributions of the first two bins, shown 
in Fig. m 

At high redshift, AGN contributions to the total in- 
frared luminosity differ from local objects. For Ltir ~ 
IO^^Lq, the AGN contribution is consistent with the lo- 
cal universe distribution given the small number statis- 
tics. For the better populated 12.5 < log(LTiR) < 13.3 
bin however, which contains most of the SMGs, the AGN 
contribution is typically lower than for the few local ob- 
jects reaching similar luminosities. All this evidence is 
tentative given the small sample size, the possible biases 
invoked by the selection and identification of the pop- 



ulations used, and the broad redshift range considered. 
It is hence too early to conclude a specific functional 
form of the change in the ^'L6agn/-^tir relation. For 
our subsequent simulations, we simply assume that the 
relation is the same as derived locally, with a change 
in the behaviour of the flat part: while in the local 
universe the increase of the AGN contribution stops at 
Ltir ^ 10^^ ® L©, for distant galaxies it occurs at lower 
luminosities (Ltir ~ 10"'^^ ^o)- It can be summarized as: 



vLq^^^I Ltir, « 



7-1.4 
''^TIR 



(Ltir < Lf) 
(Ltir > Lf) 



(1) 



Lf = 



1012-8 Lq (z<0.5) 
10^2 Lq (z>0.5) 

2.5. Discussion 



The different role of AGN at low and high redshift 
for objects in the same Ltir class should not be a sur- 
prise. Fig. [7] can be read in two (equivalent) ways: 
(1) the contribution of AGN to Ltir is lower at high- 
z with respect to the local universe for sources with 
Ltir ^ lO^^L©, or (2) the AGN content of high-z sources 
with Ltir ^ lO^^L© is similar to that of local galaxies 
with Ltir ~ IO^^Lq. This latter conclusion points out 
an analogy in the properties of SMGs and local ULIRGs. 

Such an analogy has been already discussed by 
iTacconi et al.l (|2006l ): from CO observations of a sam- 
ple of 14 SMGs, they conclude that the density of the 
molecular gas, the luminosity surface density and the 
temperature of the dust are the same in the two pop- 
ulations and SMGs are similar to local ULIRGs merg- 
ers, suitably scaled for their larger masses, luminosities 
and star formation rates, as well as their greater gas 
fractions. In spatially resolved observations, SMGs and 
ULI RGs show considera ble similarities in their kinemat- 
ics ([Tacconi et al.l l2008'): the SMGs presenting multiple 
components are interacting systems and are similar to lo- 
cal double-nucleus ULIRGs, while SMGs showing char- 
acteristics of a rotational star-forming gas disk present 
the same surface and volume densities, for example, of 
the extremely compact ULIRG Arp220. All SMGs stud- 
ied with sub-arcsecond millimeter interferometer so far, 
appear to be major mergers in different stages, similar 
to local ULIRGs. 

Our analysis, nevertheless, is model-oriented: its pur- 
pose is not to find a theoretical explanation for the dif- 
ferent populations observed, but to develop an evolution- 
ary model based on observations capable of reproducing 
available data and to help plan forthcoming surveys with 
new instruments. 

3. THE MODEL: NUMBER COUNTS AND REDSHIFT 
DISTRIBUTIONS OF INFRARED SOURCES 

3.1. The strong evolution of infrared galaxies: 
observational evidence 

There has been, in the past years, strong observational 
evidence indicating high rates of evolution for infrared 
galaxies. 

First, galaxy evolution can be observed through its im- 
print on the far-IR extragalactic background. Weakly 
constrained even as recently as the end of the 1990s, var- 
ious observations now measure or give upper/lower lim- 
its on the background from the ultraviolet (UV) to the 
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millimeter waveband (e.g. iHauser fc Dwek Il200ll ). Data 
show the existence of a minimum between 3 and 10 /Ltm 
separating direct stellar radiation from the infrared part 
due to radiation re-emitted by dust. This re-emitted 
dust radiation contains a comparable integrated power 
as the optical/near-IR: this amount is much larger than 
what is measured locally (~ 30 per cent). The CIB is 
thus likely to be dominated by a population of strongly 
evolving redshifted infrared galaxies. Since the long- 
wavelength spectrum of the background is significantly 
flatter than the spectrum of local star-forming galaxies, 
it constrains the far- IR radiation production rate history 
(|Gispert et al.|[2000f) . The energy density must increase 
by a factor larger than 10 between the present time and 
redshift z ^ 1 — 2 and then stay rather constant at higher 
redshift (tiU z - 3). 

Secondly, numerous deep cosmological surveys at 15, 
24, 70, 90, 170, 850 and 1300 /zm have resolved a frac- 
tion of the CIB into discrete sources. For all surveys, 
number counts indicate a very strong cosmological evo- 
lution of infrared galaxies, not only in the total power 
radiated but also in the shape of the LF. This is par- 
ticularly obvious at submillimeter wavelengths, where 
the background is dominated by high-luminosity galaxies 
(SCUBA and MAMBO sources). The high rates of evolu- 
tion exceed those measured in other wavelength regimes 
as well as those observed for quasars and active galactic 
nuclei (AGN). 

Finally, high rates of evolution are suggested 
by the detection of Poissonian fluctuations of the 
CIB at a high level at 60 and 100 /im with 
IRAS (iMiville-Deschenes et a l. 2002 ), and 170 /xm wit h 
ISOPHOr lfLagache et aLll2000. : MatsuIiajaj3TaI|[200l 
For ex ample, the constraints given by iMatsuhara et all 
(|2000f ) on the galaxy number counts indicate the exis- 
tence of a strong evolution in the counts. 

3.2. SEDs 

Most existing backward evolution models (§[!]) sepa- 
rate the infrared sourc es into dif ferent populations. For 
example, Franceschini at al.l (jl98 8) consider three differ- 
ent populations: (1) normal late-type galaxies, (2) in- 
teracting/starburst galaxies and (3) galaxies with AGN. 
IRAS studies showed that galaxies of different nature in 
the local Universe have different spectral energy distri- 
butions (SEDs), e.g different values of the /60//25 ra- 
tio, so usually each population is associated with a par- 
ticular SED family. The different pop ulations of the 
model ca n evolve as a single population (|Xu et al.lll998t 
IXull200df) or at di fferent rates, as suming a local LF for 
each component (jXu et al.l 1200 ll ). Nevertheless, it is 
known that infrared galaxies, in particular ULIRGs, can 
host both starburst act ivity and an AGN (|Genzel et al.l 
Il998t iLutz et al] Il998l ) and that galaxies classified as 
Seyfert show sta. rburst tracers at infrared wavelengths 
(|Lutz et al.l 12004 . while most of the f ar-IR emission m 
local QSOs is powered by starbursts (jSchweitzer et al.l 
l200fih . 

The SEDs can either be constructed empiri- 
cally, considering the individual components respon- 
sible for__tlie__overallj3mission from the galaxies 
(e.g. iRowa n- Robinson 1992), or using radiative transfer 
models fe.g. ,Efstathiou et al..,2000, ). While a complicated 
array of dust properties contributes to the SED of each 



galaxy, studies of IRAS galaxies have typically reduced 
the description to a best-fit single dust temperature, Td, 
with a one-to-one mapping to the /eo/Zioo flux ratio 
(hereafter referred to as i?(60, 100)). It has been demon- 
strated, as well, that local IRAS galaxies exhib it cor- 
relati ons of i?(60, 100) with luminosity (e.g. Dale et al.l 
I2OOID and that i?(60, 100) changes systematically over a 
large range of luminosities. A statistical relation exists 
between i?(60, 100) and infrared luminosity, even if the 
distribution is broad. Up to date, studies of the dis- 
tribution spread in temperature for each luminosity bin 
exist, but analyses of the evolution are only preliminary, 
due to shortage of information o n the temperature of 
the dust of high re dshift sources (jChapman et al.|[2003t 
IChapin et al.l l2009l ) . Difficulties of models with simple 
SED families and evo lution functions be came obvious 
with Spitzer data, (e.g. iLagache et aLll2004D . There is al- 
ready evidence for a change in SED properties with red- 
shift, with the intrinsic SED shapes of high z SMGs re- 
sembling lower luminosity local objects rather than simi- 
larly luminous local HYLIRGs, but a deep study of SEDs 
and temperatures of high- redshift galaxies is still missing. 
The main problem s are the degener acy of dust tempera- 
ture with redshift (jBlain et al.ll200^ , the limited sample 
of sub- millimeter galaxies with measurements of spcctro- 
scopic redshifts (e.g. Ch apman et al. 2005; Valiante et £d] 
l2007HPope et al.|[2008l ). and the selection effects on dust 
temperature induced by current methods to select high 
redshift infrared galaxies. 

In building the SEDs for our model, we will take into 
account the coexistence of starburst and AGN in the 
same galaxies, their varying contribution with Ltir, and 
the spread in the luminosity-temperature (L — T) rela- 
tion. AGN are hosted by infrared galaxies of different 
luminosities and their contribution to the total infrared 
emission can be either negligible or predominant. We 
studied and discussed the problem of the AGN contribu- 
tion in §[21 We found that, on average, the contribution 
to Ltir due to an AGN in IRAS galaxies is proportional 
to of the host. Even if the dispersion in the relation 
allows the presence of low luminosity AGN-dominated 
objects (e.g. NGC1068, \og{L/LQ) = 11.3), as weh as 
ULIRGs powered by pure starbursts (e.g. IRASF19297- 
0406, log(L) = 12.4 L0), this correlation holds over two- 
orders-of-magnitude, log(L) ~ 11 — 13 Lq. In § [51 we 
estimated the AGN contribution measuring the emission 
at 6 /im, after the subtraction of a starburst-like spec- 
trum; this method does not require to spatially resolve 
the AGN within the host. 

In the simulation, we go in the opposite direction. 
Each source is assigned a redshift and a total infrared 
luminosity, Ltir, consisting of both the starburst (.Lsb) 
and AGN (Lagn) contributions (Ltir = -^^sb + -^^agn)- 
Starting from Ltir, we calculate the expected value 
(i^Lguiii/LTiR) for that luminosity, using Eq. [1] (see also 
Fig. [2). 

Then, we generate a random value of i^ie^imZ-^TiR 
assuming its distribution is a gaussian with mean 
(zyLgpni/LTiR.) and sigma corresponding to the luminos- 
ity bin where Ztir lies (see Tab. [T]). From this random 
value we derive /eAGN for the redshift of the source gener- 
ated. Last, we scale an AGN template in order to achieve 
the desired fiux density /gagn at 6 fim. The AGN tem- 
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plate used is the model calculated by lEfstathiou et al.l 
(|1995f ) for the nuclear infrared continuum spectrum of 
the Seyfert galaxy NGC1068. The luminosity of the 
scaled template, Lagn, is the contribution of the AGN 
to Ltir- 

Spectral t empla tes for st arbursts are taken fro m the 
iDale et al.l (|2001l ) (see also iDale fc Heloul 120021 ) cata- 
log, divided into 64 classes from i?(60, 100) = 0.29 to 
1.64, corresponding roughly to single-component dust- 
temperature models of 23 — 45 K. SEDs are normal- 
ized by integrating each spectrum over the 8 — 1000 /im 
range and scaling to the expected value of Ltir for 
its .^(60, 100 ), calculated from the relation found by 
IChapin et al.l |2009) and adapted to the TIR luminos- 
ity (see Fig. [5]): 

i?(60, 100) - 10'='- X (^1 + X (^1 + ^ (2) 

with 7 = 0.21, S = -0.12, = -0.50 and L' = 5.6 x 
10^ Lq. In order to fit this model we follow an identical 
methodology to 'Chapin et al.' (2009"), calculating Xtir 
by integrating the same Dale & Hclou (2002) SED cat- 
alogue, normalized to IRAS 60 and 100 /im data. In this 
way we obtain an SED family in the range 10^ — 10^^'^ Hz 
spanning luminosities log(L) = 8 — 13.5 L©. During the 
simulation, in Eq. |5| we use the starburst contribution 
LsB to the total infrared luminosity only, which we ob- 
tain by subtraction Lagn from Ltir. From Eq. |2| we 
obtain the corresponding mean (i?(60, 100)) for the Lsb 
of the generated source. Then, we calculate the spread 
in the L — T relation and generate a random value of 
Rj QQ, 100), following th e gaussian distribution calculated 
bv lChapin et al.l (|2009f ) (see their Eq. 8 and Eq. 10) and 
centred on (i?(60, 100)) (see Fig. jS]). Finally, we rescale 
to the correct Lsb the lDale fc Helou I (|2002D library SED 
associated to the i?(60, 100), which was obtained from 
the gaussian distribution implementing the spread in the 
L — T relation. 

To obtain the final SED, we add the two spectra de- 
rived for the AGN and starbursts parts. Fig. |9| shows 
examples of SEDs with different Ltir and different AGN 
contributions. Because the model takes into account 
both the average relation and its dispersion, it permits a 
wide range of combinations of Ltir and AGN emission. 

3.3. Model parameters 

In this section we describe the algorithm to coherently 
model the number counts in different bands. A backward 
evolutionary model requires some basic "ingredients" : 

• the spectral energy distributions, already described 
in § 1121 depending on Ltir, i?(60, 100) and AGN 
content 

• a cosmological model describing the accessible vol- 
ume in a survey solid angle as a function of redshift. 
This is now well co nstrained from WMAP results 
(jSpergel et al.ll2007l ). A concordance cosmology of 

Ho = 75kms"^Mpc"\l^M = 0.i,flA = 0.7 is as- 
sumed. 

• a luminosity function at z = 0, e.g. the most recent 
derivation of infrared lurn i nosity function for IRAS 
galaxies bv lChapin et al.l ()2009D 



• the evolution functions for density and luminosity 
evolution. 

The solid angle observed (and therefore the number 
density of the sources) is geometrically defined by the as- 
sumed Universe model. The predicted number of sources 
in a given redshift interval (z — 0.5 dz, z + 0.5 dz) and in a 
given infrared luminosity interval (L — 0.5 dL, L -|- 0.5 dL) 
is given by 

dN{L,z)^g{z)cj,{L/fiz))^dLdz (3) 

where (j) is the local luminosity function and f{z) and g{z) 
are respectively the luminosity evolution function and 
the density evolution function and dV/ dz is the comoving 
volume element. 

The luminosity function assumed at z = is a power 
law parametr i zation of IRAS galaxies. As explained by 
IChapin et al.l ()2009[ ). the LF calculated accounts for a 
temperature bias in the underlying 60 ^m flux-limited 
sample, as well as luminosity evol ution. The LF was 
computed using the 1/Kiax method (|Schmidtl ll968') with 
a modified formalism. The parametric form of the LF for 
the Ltir is given by: 

<P{L) = ln(10)Lp, X + £) ^ 

with a = 2.42, /3 = 2.61, = 6.6 x lO^^^ Mpc^^ 
and = 6.9 x 10^° Lq (see Fig. [TUl). As in Eq. Hj we 
have re-calcu l ated t he fit for Ltir, since the results in 
IChapin et all ([200I are for narrower-bandwidth FIR lu- 
minosities. 

There exist observational evidences for strong evolu- 
tion in the infrared galaxy population (see i; 13. ip . The 
amount of observed star formation shows that a large 
component of the infrared population has evolved to the 
present epoch. Possible scenarios are density evolution, 
where, due to merging/interactions, galaxies were more 
numerous in the past, or luminosity evolution, where, due 
to enhanced star formation, galaxies were more luminous 
in the past. These two schemes of evolution capture in 
a way that is practical for computational aspects what 
in principle is a more general evolution of the luminos- 
ity function with redshift. In general, some evolution is 
needed for all populations of active sources (e.g. QSOs, 
radio galaxies and starburst galaxies) in order to model 
extragalactic source counts successfully. One of the ways 
to model luminosity and/or density evolution is using 
a simple power law. This assumption is motivated by 
similar e volutionary mode ls at radio and X- r ay wave- 
lengt hs (|Bovle et al.l Il988t iBenn et al.l 119931 : ICondon I 
Il994f ). Similar evolution has also been observed in op- 
tically selected starburst galaxies (jLillv et al.lll996f ) and 
in the submillimeter emission f rom radio loud galaxi es 
observed at 850 /xm by SCUBA (jArchibald at aliboOlD . 

Our model assumes power law luminosity, f{z), and 
density, g{z), evolutions. The following functional forms 
are adopted: 

■'^''> ~ I (1 + z)"^ (z > zi) 
5^ = 1 (1 + zr^ (z>Z2) 
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Both functions are assumed to be continuous at zi and 

The algorithm calculates the number of sources pre- 
dicted in each redshift bin (see Eq. [3]), introducing a 
poissonian scattering to the expected value. Then, a 
value Ltir is associated to each source, by Monte-Carlo 
extraction, according to the luminosity function calcu- 
lated for the redshift of the source, (j){L/f{z)). Ltir 
is then split into two parts: the fraction due to AGN 
contribution (Xagn), calculated following Eq.[l] and the 
remaining part associated to starbursts (^sb)- The SED 
is thus "built" (see § 13. 2[) adding a scaled AGN template 
(NGC1068) of luminosity Zagn with the starburst tem- 
plate corresponding to ige, after the application of a 
scatter in the L — T relation, and scaled to the appro- 
priate luminosity distance. The sources' flux densities 
in different bands are then calculated by convolving the 
redshifted SED with the bandpasses of the filters. In 
this way, we effectively simulate a virtual sky for a given 
evolution model. For each source, we know not only the 
flux densities for several bands and the ixiRj but also 
the temperature of the dust and the ise/^AGN ratio. 
The latter information, in particular, allows us to follow 
the co-evolution of starbursts and accretion, exploring 
photometric indicators to select different populations of 
infrared sources. 

4. RESULTS AND COMPARISONS WITH AVAILABLE 
SURVEYS 

Backward evolution models for the infrared spectral 
range include a significant number of tunable parame- 
ters, and using the proper observational constraints to 
adjust these parameters is the key to the successful use of 
such models. Traditionally, the first quantities to fit are 
the total IR/submm background as measured by COBE 
(GIB), and the number counts at mid-IR to submillime- 
ter wavelengths. With improving identification and fol- 
low up of SGUBA, ISO and Spitzer sources, redshift dis- 
tributions of different IR-selected populations now pro- 
vide further powerful constraints. We are now about to 
enter the next level of physical characterization includ- 
ing determination of full rest frame far-IR SEDs, a field 
that will be given a big boost with BLAST and Herschel, 
and better disentangling the role of star formation and 
AGN in high redshift infrared galaxies. The model pre- 
sented here is designed to allow comparison with this new 
type of constraints on SED and AGN content. By implic- 
itly coupling AGN evolution to infrared galaxy evolution, 
this model is also amenable to comparison with AGN sur- 
veys at other wavelengths, if making use of assumptions 
about AGN SEDs and obscuration. 

The most relevant results of the last years are related 
to surveys carried out by SCUBA and Spitzer. Thus, 
we will compare our results with the data available for 
these instruments . H e reafte r , we will use number c ounts 
from IPapovich eFn] (I2004D and 'Shupe et all (pOOS ) for 
the 24 /n m sources, " iFraver et al] (j2006a .lil for the 70 
sources, IFraver et al.l (|2006al ) for the 160 /im sources, 
ICoppin et al.l (|2006D for the 850 ^m s o urces . Red- 
shift distributions are from ^ uyts et all (|2008D (24 /im 
sources) and Cha pman et al.l (j2005) (850 /im sources). 
Some of the models si mulated are a l so compared with the 
GIB measurements ()Fixsen et al.l 119981 : iLagache et al.l 



20001 iRenault et aL l lMll : iMiville-Deschenes et all 120021 : 
Lagache et al. All the simulations are run on 

a field of view of lOdeg , in order to be comparable 
with the main infrared surveys and sufficiently sample 
the high luminosity end. The redshift range covered is 
< z < 8. In the plots of the redshift distributions 
(see Fig. [TT]-Fig. [H]) . the available data have been ap- 
propriately scaled in order to be comparable with the 
simulations. 

In order to find the best set of parameters 
(rii, 712, zi,mi, 7712, Z2) fitting available data we have run 
simulations on small fields of view, changing slightly each 
parameter "by hand" . Changing the parameters one by 
one allows us to understand how each parameter influ- 
ences the final results: we have used this approach to 
have full control of the simulation and to understand 
the role and weight of each parameter. During the opti- 
mization process, we have considered differently the sev- 
eral constraints given from current surveys, depending on 
their reliability: first, we have worked to reproduce the 
most accurate measurements and only later we have tried 
to fit the less constra ined results. For example, SHADES 
(jCoppin et al.l l2006l ) has made accurate measurements 
of the number counts at 850 /im down to faint fluxes 
(~ 2 mjy) and Spitzer has provided very sensitive mea- 
surements with minimum uncertainties of the number 
counts at 24 /xm, either with observations on large scales 
(49 deg2 , iShupe et al.ll200l or ve ry deep, down to fluxes 
~ 0.1 mJy (Papovi cii et al.ll2004f l. On the other hand, 
the redshift distribution of submillimeter galaxies is not 
well cons trained at high 2 : . In fact, the radio-selected 
sample of lChapman et al.l ()2005D is biased towards z <'i 
because of the radio flux limit, and there are severa l hints 
of ga l axies at z ^ 4 (Valiante et al. 2007; Knudsen et al.l 
120081 iDaddi et al.]|2009lV The redshift distribution of 
24 /im sources are obtained from photometric redshifts 
on a jvT.s-band sele cted galaxies in the GOODS-S field 
(|Wuvts et al.ll200il : cosmic variance can have an impor- 
tant role on the redshift distribution in a field so small, as 
well as the selection and the choice of the SEDs for pho- 
tometric redshifts. All these effects could overestimate 
the numbe r of low redshift obj ects. Finally, the surveys 
at 160 /xm (jFraver et al.ll200"6af l are still rather shallow. 

Below, we present a summary of the optimization pro- 
cess, showing the main steps followed to obtain the best 
parameter set (ni, n2, zi, mi, m2, Z2). A summary of the 
parameters for all the models discussed are shown in 
Tab.m 

4.1. Pure starbursts SEDs 

As a starting point, we compare the observations with 
a model considering only the starburst part of the SED, 
including its scatter in the L — T relation (model Ml). 
The submillimeter galaxies are starburst dominated, so 
we look for the best parameter set reproducing at least 
the 850 /xm number counts. We assume a strong evolu- 
tion in luminosity (rxi = 3, 7x2 = 0) and a weaker evolu- 
tion in density {mi — 1, 7712 — 0), with 2:1=22 = 2. The 
results are shown in Fig. [TT] The bright 850 /xm counts 
are well reproduced. However, all the other counts are 
underestimated: using this evolutionary model, there is 
still room to introduce an SED taking into account the 
AGN contribution. The redshift distribution at 24 /xm 
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is quite well reproduced, while at 850 /xm the modeled 
redshift distribution shows a strong excess at z > 3.5. 
Adding the AGN part, we can be confident that the 
submiUimeter distribution will not change much, but we 
should not expect the same trend in the mid-IR, since 
there the AGN contribution can be significant, in partic- 
ular at high redshifts. 



4.2. Adding the AGN contribution 

In this second step, we build SEDs as described in 
§ 13.21 In particular, we determine the AGN contribu- 
tion following Eq. [TJ As calculated in § [2l we assume a 
gaussian distribution in each infrared luminosity bin (see 
Tab.[T]). The results, assuming the same evolution ap- 
plied previously, are shown in Fig. [T^] (M2). Compared 
to model Ml, the 24 ^m number counts are now better 
reproduced, even if the shape of the peak is slightly dif- 
ferent. The redshift distribution at 24/ini has become 
worse: the AGN contribution has introduced a large ex- 
cess at z ~ 2. The parameter set used up to now is no 
more acceptable and we have to explore the parameter 
space further, in order to reproduce all the observations. 

Comparing different kinds of evolution, we learn that 
we need a strong luminosity evolution in order to repro- 
duce the number counts at 24 /im and the 850 /im, but 
this always introduces an overprediction of the 24 /im 
z ^ 2 sources. On the other side, with a strong density 
evolution we obtain a redshift distribution at 24 /xm that 
is very similar to the observed one, but we completely 
lose the peak in the 24 /im number counts. We note, 
moreover, that a shift of the peak of the density evolu- 
tion, Z2, towards lower redshifts, causes a change in the 
shape of the redshift distribution at 24 /im, reducing the 
number of the z ~ 2 sources, without changing the 24 /tm 
number counts. Last, if we want to reduce the number 
of high redshift 850 /tm sources, we need a luminosity 
evolution peaking at zi and then dropping down (n2 has 
to be negative). Fig. 1131 shows number counts and red- 
shift distributions for a model with ni — 3.4, rt2 — —1, 
zi — 2.3, mi = 1, 7712 = —1.5 and Z2 = 1 (M3). Data 
at 24 /im are well reproduced, in particular the predicted 
peak of the number counts has now the same shape of the 
measurements even if it is slightly lower, while at 850 /tm 
the simulation underpredicts the integral number counts 
but reproduces the redshift distribution. 



4.3. Evolving the L — T relation 

In this simulation, the starburst part of the SED still 
does not vary from the local to the distant universe. This 
last point is in contradiction to recent works already 
mentioned in § 13.21 showing that submiUimeter galax- 
ies, even being as luminous as local "HYLIRGs", have 
dust temperatures typical of less lumin ous local galax- 
ies (ULIRG s) (Chapman ct al. 2005; P ope et al.ll2006l : 
IValiante et all 120071) . IPope et al.l ([2001), in particu- 
lar, needed to modif y the "classical" SED model of 
iCharv fc ElbaJI ([20?)lh increasing the cold component of 
the dust, in order to fit the SEDs of the submiUimeter 
galaxies in the GOODS-N field. 

We then include in our model a simple evolution in the 
L — T relation, accounting for the cooler temperatures of 
high redshift infrared galaxies. We modify Eq. [2] in the 



following way: 
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with all the parameters as in Eq. [51 This evolution is 
consistent with the observations cited above. Fig. [14] 
shows number counts and redshift distributions for this 
new model with ui = 3.4, n2 = —1, zi = 2.3, mi = 1, 
m2 = — 1.5 and Z2 — 1 (M4). All the submiUimeter mea- 
surements are now well reproduced, with perhaps some 
excess in the number counts at strong fluxes, together 
with the mid-IR data, whose number counts have im- 
proved and redshift distribution has not felt the effect 
of the new SEDs. The number counts in the far-IR, in 
particular at 160 /im, are nevertheless underestimated. 

Even if the peak of the 24 /tm number counts is well re- 
produced by the model M4, there is still some excess in 
the number counts at strong fluxes, where AGN dom- 
inate, but the uncertainties are large given the small 
number statistics. We compare the number of AGN pre- 
dicted by our model with the number of sources expected 
adopting QSO luminosity functions. We select only the 
brightest 24 /im sources (/24 > 5 mJy) and consider only 
the AGN-dominated sources (iACN/^SB > !)■ With the 
exception of some local and faint objects, all the sources 
selected in this way are highly luminous, distant AGN. 
Applying a simple bolometric correction, Lbol ^ IQvLi, 
at mid-IR rest wavelengths, we can calculate the ex- 
pected number of QSOs for the same Lbol and z ranges 
of our subsample . We integrate the luminosity function 
of lHopkins eTall (jlOOIi) ("full" model, with pure lumi- 
nosity evolution and bright- and faint-end slopes evolv- 
ing with redshift) for 10^2.2 < j^^^^ < ;^q14.5 j^^ 

0.1 < z < 2.3. The number of the sources selected 
from our simulation does not exceed the expected num- 
ber given by the bolometric luminosity function. Never- 
theless, we already pointed out that the high part of the 
i^L6(AGN)/Ltir distribution is not well constrained (see 
§ 12. 3p . In particular, the hypothesis of the gaussian dis- 
tribution could be wrong: it was assumed following the 
principle of the minimization of the free parameters, but 
a different kind of distribution, for example an asym- 
metric distribution with a cut-off at high luminosities, 
cannot be excluded (see § [6]) . 

4.4. Evidence for a local "cold" population? 

Our models were no able to reproduce the num- 
ber counts at 160 /im, underestimating them by a fac- 
tor of ~ 5. The ISOPHOT Serendipity and FIR- 
BACK s urveys have re veale d a po pulation of nearb y cold - 
galaxies (IStickel et al.lll998ll2000l: IChapman et al.ll2002l: 



Patris et al.l 120031: iDennefeld et al.l 120051: ISaiina et al l 
2006( 1 ■ under-represented in the 60 /tm IRAS sample. 
These objects are often associated with bright optical 
spiral galaxies, and their far-IR colours (/170//100 1-3) 
indicate a rising spectrum beyond 100 /tm, similar to 
that seen for exa mple in the Milky Way galactic ridge 
(ISerra et al 11197^. 

Lagache et al.l ( 20031 1 implemented this class of objects 
(hereafter "cold galaxies") in their evolutionary model: 
cold galaxies dominate the z = luminosity function 
at low luminosities and become less important at higher 
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redshifts, because their evolution is passive and short 
(ni — 1,712 — 0, zi — 0.4). Their contribution to the 
high luminosity part of the luminosity function is still 
substantial up to £tir ~ l O^^Lfri, while at larg er lu- 
minosities it is quite small. La gache et al.l ([2003) show 
that the fraction of cold galaxies can contribute to the 
total number counts up to ~ 50% at 170 \m\ in the flux 
range where a comparison with measurements is possible 
(0.1 — l.OmJy, see their Fig. 9). The modifications added 
later to th e model have not cha nged the number counts 
at 160 /im (jLagache et al.l l2004 their Fig. 4). 

In order to introduce a similar population in our model, 
we have to assume that probably the s catter in the L — T 
relation calculated by IChapin et al.l (|2009[ ) and imple- 
mented in our simulations needs some corrections in or- 
der to take into account the cold galaxies population. 
The lack of these corrections can explain why we un- 
derestimate the 160 ^m number counts. Moreover, any 
changes of the adopted L ~ T relation and its scat- 
ter has to be constrained by the 850 /im counts. The 
160 /im counts need to be increased without creating ex- 
cess 850 /im counts. This may support a mostly local 
change. Moreover, colder galaxies with a higher 160 /xm 
flux have typically lower 70 /xm flu x: any change intro- 
duced cannot be too much strong. IChapin et al.l (120091 ) 
discussed two SMGs that are outliers respect to their 
L — T relation: they explain the exceptions assuming a 
mismatching in the optical counterparts. Nevertheless, 
they do not exclude that both the galaxies are really at 
low redshift and a population of cold galaxies do exist at 
z < 1 and are largely missed in IRAS surveys. 

We then modify the spread in the L — T relation as- 
suming an asymmetric gaussia n distribution . The orig- 
inal distribution calculated bv IChapin et al.l ()2009[ ) is a 
standard gaussian in log(i?(60, 100)). We assume a 7 
times broader gaussian on the "cold" side of the peak 
respect to the "warm" side at low redshifts (z < 1). In 
this way more cold galaxies are generated with respect 
to the IRAS population. The results for this model with 
ni — 3.4, — — Ij zi = 2.3, mi — 2, TO2 — —1.5 and 
Z2 = 1 (M5) are shown in Fig. [TS] Finally, both the 
far-IR (160 /im and 70 /im wavebands) number counts 
are reproduced within the uncertainties as well as the 
submm measurements. We note that a qualitatively sim- 
ilar need to add local cold sources arose when using in 
our models both the ear lier IRAS-based determ ination of 
temperature sp read by iCha pman et all ()2003f ) and the 
recent work of IChapin et al., (2009. ) . We can consider 
M5 our best-fit model: all the number counts are well 
reproduced (within a factor < 2), and also the redshift 
distributions are very similar to the observed ones, even 
if there are some discrepancies, due probably to the rea- 
sons explained above. 

An example of how this population of cold galaxies 
could influence results of Herschel surveys is given in § [5] 
(see Fig.^ottom). 

A further constraint on our evolutionary model comes 
from the comparison with the current CIB measure- 
ments. Our model predictions for the CIB are derived 
by summing up the flux densities of all sources for the 
bands in question. FiglTOIshows the predicted CIB inten- 
sity at specific wavelengths together with the comparison 
with present observations for the model M5: it agrees 
with all the available measurements and limits. We also 



compare the IR luminosity function calculated at dif- 
ferent redshifts with measurements. Fig. [T7] shows the 
luminosity function of our be st model M5 fitting all data 
from low redshifts z = 0.3 (iHuang et al.ll2007l) to z = 
0.4, 0.5,0.7,0.9,1.1 measur ed by [LeFloclieial. (200^ 
and iMagnelli et all (|2009D . z = 1,2 bv ICaputi et all 
(12007D and SMGs at z - 2.5 (|Chapman et al.ll2005l) . 

5. STAR FORMATION HISTORY 

There has been much interest in the star formation 
rate and met a l prod uction history of the Universe since 
iMadau et al.l (|l996l ) showed that they could be derived 
from deep UV/optical cosmological surveys. Since then 
the so-called Madau diagram has been revised many 
times through various improvements, including (1) the 
corrections for the effec t of dust exti nction on UV/optical 
luminosities fe.g. lSteidel et al.lll999( ). (2) result s from less 
extin ction-sensitive Balmer line surveys (e.g. lYan et al] 
Il999f ). an d (3) results from sub millimeter SCUBA sur- 
veys (e.g. IChapman et al.|[2005l ). 

Adopting the conversion factor of iKennicutt I ()1998D . 



SFR [Mo yr-i] - 4.5 • 10-*^ x Ltir [erg s'l] 



(8) 



we can convert the cosmic luminosity density evolution 
of our model to a SFR curve. In the calculation, we take 
into account only the fraction on infrared luminosity due 
to star formation (Lsb)- In Fig. [18] the results of the sim- 
ulation using model M5 are compared with the survey 
data already shown bv IChapman et aTl (|2005l ) (see their 
Fig. 12). The model is in very good agreement with the 
observations for the full redshift range. In particular, for 
z > 4, the model results are similar to the results from 
SCUBA surveys, while it is slightly lower than the mea- 
surements obtained from the UV/optical surveys. This 
does not surprising, since the model itself is constrained 
by the infrared/submillimeter data. 

6. AGN CONTRIBUTION TO THE INFRARED EMISSION 

As discussed earlier, LIRGs and ULIRGs, with their 
huge infrared luminosities, correspond to an extremely 
active phase of dust enshrouded star formation and/or 
AGN activity. They become an increasingly significant 
population at high redshift, representing an important 
phase in the buildup of massive galaxy bulges and in 
the growth of their central supermassive black holes. To 
understand these processes, it is essential to separate the 
relative contributions of AGN and starburst activity to 
the infrared luminosity of LIRGs and ULIRGs. 

Recent work has shown how the IRAC color-color di- 
agram and the MIPS 24 to 8 /xm color can be used 
to ident ify AGN-d ominated so urces (Lacv et al. 2004; 
Saiina et all |200^ iStern at al.l [20051: lYan et all 12004 ; 
Brand et al.[|2006[ ). iBrand et al.l (|2006f ). in particular, 
have demonstrated how the 24 to 8 /im flux ratio (C = 
log[7^/„„(24/im)/i//,y(8 /im)]) can be used to disentangle 
the contribution of AGN and starbursts to the total re- 
processed mid-IR (~ 5 — 25 /im) emission as a function 
of the 24 /im flux. We use this study as a first exam- 
ple of how to compare our model with observations that 
constrain the AGN content of high— z sources. 

Nevertheless, a photometric study of this type is sub- 
ject to several caveats. Many broad emission- and 
absorption-line features are known t o be present in the 
mid-IR spectrum of ULIRGs (e.g. iHouck et all 120051 : 
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lYan et al.|[2005D . and this may affect C as function of 
redshift. For example, at 1-1 ^ 2 < 1.7, the 24 /xm ob- 
served emission can be strongly attenuated by the sili- 
cate absorption feature at 9.7 fim. It is also possible that 
an AGN could be heavily embedded in large amounts 
of cooler dust and could remain undetected in the 8 /im 
band even though it dominates the bolometric infrared 
emission. In practi ce: on one hand, the galaxies, which in 
iBrand et aT] ()2006[ ) are classified as "AGN-dominated" , 
do host an AGN, but it is not necessarily the dominant 
source powering the infrared emission; on the other hand, 
objects that are known to be dominated by AGN in the 
mid-IR, like NGC1068, do not present the color (C ~ 0) 
that i s considered " typica l" for AGN-dominated sources 
in the IBrand eTell (|2006[ ) study. 

Obtaining mid-IR spectroscopy of complete unbiased 
24 /im selected samples of infrared sources will be impor- 
tant in putting more constraints on the AGN contribu- 
tion to the total infrared emission and in characteriza- 
tion of the 24 fim selected population. Large Spitzer-IRS 
observing programs covering sources with 24 /im fluxes 
around and below 1 m Jy are under analysis and will be 
well suited for this task. In combination with IRS sur- 
veys at brighter fluxes, they will trace the evolution of 
the AGN/SB ratio, strength of PAH emission and mid-IR 
opacities as a function of L^m and z. In t he meantime, 
result s of photometric studies like the one of lBrand et al.l 
(l200l can be compared to our models (see Fig. [!%,). 

The models developed here allow us to predict the frac- 
tion of AGN-dominated sources as a function of the mid- 
IR flux. In particular, it is possible to distinguish be- 
tween sources dominated by AGN at 24 /zm and sources 
AGN-dominated with respect to the total infrared lumi- 
nosity. This distinction is important when comparing the 
results of the simulation with the observations: the infor- 
mation about /24 will be useful in comparison with the 
upcoming results from the IRS observations cited above, 
while the predictions related to Xtir will be confronted 
with Herschel results. 

At the moment, as already explained, the shape of the 
AGN distribution is not well constrained at the highest 
ratios of AGN to total infrared luminosity. As an exper- 
iment, in order to reduce the number of luminous AGN, 
we change the assumed gaussian distribution introducing 
a cut for large values of i^ie/im/^TiR: we reject all val- 
ues of J^ie^m/^TiR more than 3 a away from its expected 
value. This model, again using the same parameter set 
(ni = 3.4, n2 = —1, zi = 2.3, mi = 1, 1712 = —1.5, 
Z2 — 1, model M6), does not significantly change the 
number counts and redshift distributions with respect to 
the best model M5, so we do not show them here. 

Differences between the two models are evident when 
doing different kinds of analysis. Fig. [19b shows the 
fraction of all sources whose mid-IR emission is dom- 
inated by AGN (/24AGN//24SB > 1) as a function of 
the 24 /im flux density for model M5 {solid histogram) 
and M6 {dashed histogram). The trend an d the values 
obtain ed are similar to those measured by IBrand et al.l 
()2006[ ) {open squares), even if the measurements have to 
be compared carefully for the reasons explained above. 
Our predictions will be easily compared instead with the 
upcoming IRS results, that will help us to put more con- 
straints on the modeled AGN contribution. In fact, the 
AGN fraction of model M6, where the AGN contribu- 



tion has been slightly modified, is a bit lower compared 
to model M5: only more constraints coming from mid-IR 
spectroscopic observations can discriminate between the 
two models M5 and M6. Fig. [19b shows the fraction of 
all sources whose total infrared emission is dominated by 
AGN (Xagn/^sb > 1) as a function of the 24 /xm flux 
density for the two models M5 and M6. As expected, the 
values of the fraction are lower, because the AGN starts 
to dominate the mid-IR part of the SED before being 
dominant in the entire infrared range (see e.g. Fig. (HI). 

7. OPEN ISSUES 

Our best-fit model M5 is able to reproduce most of 
the measurements available. Nevertheless, there are still 
some issues that need a deeper analysis and discussion 
and will be dealt with in future studies. 

The first point concerns the so called "local cold pop- 
ulation" . Even if there are several clues about the exis- 
tence of a population of colder IR galaxies not detected 
by IRAS, only forthcoming surveys at longer wavelengths 
carried out with Herschel will be able to reveal and char- 
acterize this kind of objects. 

The second point regards the distribution of the AGN 
contribution and the fraction of AGN-dominated sources. 
We already pointed out that the models M5 and M6 pre- 
dict similar 24/im number counts, but different amount of 
AGN-dominated sources. Galaxies at z ~ 2 presenting a 
"mid-IR excess" have redder K—5.8 /im colours than nor- 
mal galaxies: this is evidence for an AGN contribution to 
the mid-IR continuum due to warm dust. The presence of 
Compton-thick AGN is confirmed by the stacked Chan- 
dra X-ray data (iDaddi et al.ll2007h . The sky density and 
the volume density of this population of obscured AGN 
agree reasonably well with those predicted by the back- 
ground synthesis models of iGilli et al.l ()2007[ ). In order 
to discriminate between our best models, M5 and M6, we 
need to compare the AGN-dominated sources predicted 
with upcoming results with IRS. This will introduce a 
further constraint, in addition with the bounds already 
used to the total number counts, and it will enable the 
final model to reproduce the co-evolution of AGN and 
starbursts in more detail. 

The third point regards the evolutions of the AGN— 
and Color— Ltir relations. It is clear that some evolu- 
tion is needed in order to reproduce existing observations, 
nevertheless the evolution implemented in this model is 
quite simple. More information regarding the AGN con- 
tent of high— z far-IR selected sources and the shape of 
the SED of SMGs will enable a better characterization of 
these details. Understanding evolution is intimately re- 
lated to the black hole growth and to the star-formation 
rate history of massive galaxies. 

8. PREDICTIONS FOR MULTIBAND HERSCHEL SURVEYS 

In this section, we show predictions using the best-fit 
model M5. We concentrate on future surveys with Her- 
schel, which will be launched in early 2009, in particular 
on the PEP survey^. All three PACS bands (70, 100, 
160 /im) are considered. Predicted number counts at 70 
and 160 /im are already shown in Fig. [14] and compared 
with the currently available data in those wavebands. 

In Fig. \20t ov the predicted integral number counts at 
100, 250, 350, 500 /im are plotted. While at 70 /im the 

^ http:/ /www.mpe.mpg.dc/ir/Research/PEP/ 
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AGN-dominated sources can give a non-negligible con- 
tribution to the total number counts, in particular for 
fro ^ 50 mjy, in the other Herschel bands the number 
counts are exclusively dominated by the starburst com- 
ponent. 

Fie. \20]m iddle shows predicted redshift distributions for 
the deep (/loo > 1-5 mJy) and large (/leo > 8mJy) PEP 
fields. The deep (10' x 15') and the large (85' x 85') 
fields are centred respectively on the GOODS-S and on 
the COSMOS fields, in order to have the maximum avail- 
ability of multi- wavelength data and follow up opportuni- 
ties. The chosen fields are in fact fully covered in several 
of the following: deep X-ray surveys, UV/optical/near- 
IR/IRAC imaging, HST imaging, Spitzer mid-IR sur- 
veys, submillimeter surveys and radio mapping. The pre- 
dicted distributions show a prominent peak at z ~ 1 and 
a shoulder reaching z ~ 2 and beyond. Both the surveys 
will observe sources up to z ~ 3, allowing a consistent 
step forward in the study of the cosmic evolution of dusty 
star formation and of the infrared luminosity function, as 
well as in the determination of the overall SEDs of active 
galaxies. 

In Fis\2Ubottom the predicted redshift and dust tem- 
perature distributions of galaxies selected in the PACS 
bands (100 and 160 /im) for models M4 (no local cold 
population) and M5 (with local cold population) are 
plotted. The selection is assumed to be the same as 
expected in the CDFS field (~ 8deg^ /loo > 27mJy, 
/i60 > 38.8 mJy). Redshift distributions are similar in 
shape, but M5 presents an increase in the number of 
detected objects at z < 1 of a factor ^ 3. Dust temper- 
ature distributions are completely different: M4 predicts 
a mean temperature of Td ~ 38 K, while M5 presents a 
skewed distribution with a prominent peak at ~ 23 K 
and a tail up to 60K. Multiband Herschel surveys 
will measure these distributions. 

9. CONCLUSIONS 

The aim of this work has been to find a "backward 
evolution" model that can viably fit the galaxy source 
counts and redshift distributions from the mid-infrared 
to submillimeter wavelengths whilst not violating the 
constraints set by the Cosmic Infrared Background mea- 
surements. Increasingly detailed observations are char- 
acterizing high redshift galaxies simultaneously at many 
infrared wavelengths, and are starting to unravel the 
role of AGN within them, suggesting a refinement of 
previous models. We adopted a Monte-Carlo-based ap- 
proach that considers the evolution of a coherent popu- 
lation of infrared galaxies with distribution functions de- 
scribing the contribution of AGN and the spread around 
the luminosity-temperature relation for the far-infrared 
emission. 

As a necessary ingredient of this model, we have char- 
acterized the local distribution of AGN contributions as 
a function of total infrared luminosity. We have applied 
spectral decomposition to a large sample of Spitzer-IRS 



spectra of ULIRGs and LIRGs with Ltir > 10^^ Lq 
to isolate the AGN 6/im continua. The distribution of 
ieAGN/^^TiR changes with Ltir. We compare the AGN 
detections and limits to simple Monte-Carlo simulations, 
making the assumption of an intrinsically gaussian dis- 
tribution, to quantify this increase with luminosity of the 
AGN contribution to the infrared luminosity. The best 
fit, vLq ^ni/ Ltiji oc Lxiri gives a — 1.4 ± 0.6. The rela- 
tion does not hold any more at high redshifts. The most 
luminous high redshift infrared galaxies, Ltir IO^^^Oj 
show a small contribution from AGN, being mainly star- 
burst powered. 

The local start of our backward evolution models 
combines this result with the lu minosity function an d 
temperature-luminosity relation of Chapin et al. (2009). 
We have explored several models considering evolution- 
ary variation of increasing subsets of the parameters 
within this framework. We find that we can satisfactorily 
reproduce number counts and redshift constraints from 
the literature with such a population of infrared galaxies 
evolving in luminosity and density. The contribution of 
starbursts and AGN varies with redshift and luminosity, 
as does the far-infrared dust temperature. The luminos- 
ity evolution is mainly constrained by number counts at 
24 and 850 /im, while the density evolution is strongly 
influenced by the 24 ^m redshift distribution. While this 
evolution in luminosity, density, AGN contribution and 
dust temperature is clearly suggested by this sequence 
of models, the detailed functional forms and parameters 
are necessarily uncertain with current data^ 

As in previous work (e.g. Lagachc ct al. 20031 ). we have 
realized a need to invoke a population of local cold galax- 
ies in order to reproduce the number counts at 160 /im. 
Even with such clues about the existence of this popu- 
lation, its observational characterization remains incom- 
plete. 

The coherent consideration of both AGN and starburst 
in our model allows for direct comparison to observations 
constraining the role of AGN in high redshift infrared 
galaxies. Trends with modelled mid-infrared flux from 
our best fit model are con sistent with photometric es- 
timates ( Brand et al.ll200"6l ). Forthcoming IRS observa- 
tions of large flux-limited samples will both reduce the 
observational uncertainties in characterizing the role of 
AGN and allow to better constrain the role of the most 
luminous AGN in our models. 

Characterization of the SEDs of z ~ 0.5 — 3 galax- 
ies near their rest frame far-infrared peak will make a 
big step forward with the Herschel mission, and advance 
the understanding of the AGN-galaxy co-evolution. The 
model type presented here is prepared for the substan- 
tial tests and refinements of assumptions on the evolu- 
tion of role of AGN and luminosity-temperature relation 
that will be possible. As a first step, we have presented 
predictions from our current best-fit model for redshift 
distributions and the role of AGN in deep surveys with 
Herschel-PACS. 
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TABLE 1 

Full sample of infrared galaxies: results of the fit and 
corrected /eagn fluxes 



Name 


2 


log{LTIR) 


/sagn 


/6AGN 








mJy 


mJy 



Sources from RBGS catalog (| Sanders et al.|[200^ ') 



NGC0023 


0.015 


11. 


05 


13.9±1.0 


<23.8 


NGC0034 


0.020 


11. 


44 


48.8±3.1 


<62.9 


MCG-02-01-051/2 


0.027 


11. 


41 


6.8±0.8 


<19.9 


ESO350-IG038 


0.021 


11. 


22 


57.8±0.5 


53.4±4.0 


NGC0232 


0.020 


11. 


30 


20.5±1.5 


<31.4 


MCG+12-02-001 


0.016 


11. 


44 


20.5±1.7 


<48.4 


NGC0317B 


0.018 


11. 


11 


14.6±0.7 


<23.6 


IC1623A/B 


0.020 


11. 


65 


196.4±2.5 


162.0±30.5 


MCG-03-04-014 


0.035 


11. 


63 


8.3±0.7 


<24.1 


ESO244-G012 


0.023 


11. 


39 


12.6±1.3 


<34.3 


CGCG436-030 


0.031 


11. 


63 


35.8±0.6 


26.5±8.2 


ESO353-G020 


0.016 


11. 


00 


20.2±1.1 


<30.1 


ESO297-G011/012 


0.017 


11. 


09 


16. Oil. 1 


<31.1 


IRASF01364-1042 


0.048 


11. 


76 


4.0±0.2 


<5.2 


IIIZw035 


0.028 


11. 


56 


2.1±0.1 


<2.7 


NGC0695 


0.033 


11. 


63 


5.5±0.5 


<13.5 


UGC01385 


0.019 


10. 


99 


7.1±0.7 


<16.5 


NGC0828 


0.018 


11. 


31 


17.3±1.1 


<29.2 


NGC0838 


0.013 


11. 


00 


5.6±1.5 


<36.0 


IC0214 


0.030 


11, 


.37 


3.9±0.5 


<11.6 


NGC0877 


0.013 


11, 


.04 


3.7±0.1 


3.3±0.4 


UGC01845 


0.016 


11, 


.07 


25.2±2.0 


<42.5 


NGC0958 


0.020 


11, 


.17 


5.0±0.1 


4.4±0.6 


NGC0992 


0.014 


11, 


.02 


5.0±1.0 


<27.6 


NGC1068 


0.003 


11, 


.27 


10493. 6±30. 6 


10482.0±32.3 


UGC02238 


0.021 


11, 


.26 


16.4±1.4 


<35.9 


IRASF02437+2122 


0.023 


11, 


.11 


18.3±0.4 


15.8±2.3 


UGC02369 


0.031 


11, 


.60 


9.3±0.6 


<15.3 


UGC02608 


0.023 


11, 


.35 


67.0±1.1 


60.0±6.3 


NGC1275 


0.018 


11, 


.20 


120.2±0.8 


119.7±0.9 


IRASF03217+4022 


0.023 


11, 


.28 


12.9±0.9 


<18.9 


NGC1365 


0.005 


11, 


.00 


227.8±1.7 


213.9±12.4 


IRASF03359+1523 


0.035 


11, 


.47 


0.4±0.1 


<0.7 


CGCG465-012 


0.022 


11, 


.15 


6.3±0.5 


<12.2 


IRAS03582+6012 


0.030 


11, 


.37 


189.7±1.5 


188.7±1.7 


UGC02982 


0.017 


11, 


.13 


6.0±0.8 


<19.4 


ESO420-G013 


0.012 


11, 


.02 


86.6±1.1 


69.5±15.2 


NGC1572 


0.020 


11, 


.24 


18.6±0.6 


<24.9 


IRAS04271+3849 


0.019 


11, 


.06 


9.3±0.7 


<21.4 


NGC1614 


0.016 


11, 


.60 


24.8±2.9 


<56.7 


UGC03094 


0.025 


11, 


.35 


21.2±0.5 


<27.2 


ESO203-IG001 


0.053 


11, 


.79 


12.9±0.4 


12.6±0.5 


MCG-05-12-006 


0.019 


11, 


.12 


13.9±0.5 


<20.3 


NGC1797 


0.015 


11, 


.00 


12.8±0.9 


<24.8 


CGCG468-002 


0.017 


11, 


.10 


41.1±0.6 


33.8±6.5 


VIIZw031 


0.054 


11, 


.94 


8.5±0.7 


<21.0 


IRAS05083+2441 


0.023 


11, 


.21 


5.1±0.9 


<28.2 


IRAS05129+5128 


0.027 


11, 


.36 


11.4±0.4 


<17.5 


IRASF05187-1017 


0.028 


11, 


.23 


6.2±0.3 


<8.5 


IRASF05189-2524 


0.043 


12, 


.11 


204.8±0.9 


201.3±3.2 


IRAS05223+1908 


0.030 


11, 


.59 


367.8±1.1 


366.7±1.5 


NGC1961 


0.013 


11, 


.02 


8.6±0.3 


7.8±0.8 


MCG+08-11-002 


0.019 


11, 


.41 


26.3±1.1 


<38.6 


UGC03351 


0.015 


11, 


.22 


22.3±0.9 


<31.8 


IRAS05442+1732 


0.019 


11, 


.25 


13. Oil. 5 


<33.2 


UGC03410 


0.013 


11, 


.04 


5.2±0.5 


<11.3 


IRASF06076-2139 


0.037 


11, 


.59 


9.7±0.3 


7.1±2.3 


NGC2146 


0.003 


11, 


.07 


45.6±6.2 


<150.3 


ESO255-IG007 


0.039 


11, 


.84 


11.7±1.3 


<36.6 


ESO557-G002 


0.021 


11, 


.19 


5.3±0.3 


<9.0 


UGC3608 


0.022 


11, 


.30 


8.0±0.4 


<12.4 


IRASF06592-6313 


0.023 


11, 


.17 


11.7±0.5 


<16.5 


AM0702-601 


0.031 


11, 


.58 


100.2±0.4 


98.0±2.0 


NGC2342 


0.018 


11, 


.25 


5.6±0.4 


<11.5 


NGC2369 


0.011 


11, 


.10 


30.0±1.4 


<42.5 


IRAS07251-0248 


0.088 


12, 


.32 


16.4±1.1 


15.9±1.2 


NGC2388 


0.014 


11, 


.23 


22. Oil. 1 


<34.5 


MCG+02-20-003 


0.016 


11, 


.08 


53.4±1.7 


46.2±6.6 


IRASF08339+6517 


0.019 


11, 


.05 


2.9±0.8 


<15.7 


NGC2623 


0.018 


11, 


.54 


22.3±0.8 


<28.8 


IRAS08355-4944 


0.026 


11, 


.56 


78.7±0.6 


69.1±8.5 
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TABLE 1 

continued - Full sample of infrared galaxies: results of the 

FIT and corrected /sAGN FLUXES 





z 


log{LTIR) 


/eAGN 




/eAGN 










mjy 




mJy 


ESO432-IG006 


0. 


.016 


11 


.02 


18. 


.7±1. 


.0 


<28.8 


ESO60-IG016 


0. 


.046 


11 


.76 


56. 


.4±0. 


.6 


53.0±3.1 


IRASF08572+3915 


0. 


.058 


12, 


.10 


291. 


.5±2. 


.9 


290.5±3.0 


IRAS09022-3615 


0. 


.060 


12 


.26 


83. 


.2±0, 


.5 


75.9±6.5 


IRASF09111-1007 


0. 


.054 


12 


.00 


5, 


.8±0, 


.3 


<8.4 


UGC04881 


0. 


040 


11 


.69 


9, 


,0±0, 


.3 


<13.6 


UGC05101 


0, 


039 


11 


.95 


46. 


7±1. 


.0 


40.3±5.8 


MCG+08-18-013 


0, 


,026 


11 


.28 


0, 


.3±0, 


.1 


<0.3 


IC0563/4 


0, 


,020 


11 


.19 


6, 


,0±0, 


.4 


<10.0 


NGC3110 


0, 


,017 


11 


.31 


9. 


,1±0, 


.7 


<21.0 


IC2545 


0, 


,034 


11 


.73 


311. 


7±6. 


.2 


308.8±6.7 


IRASF10173+0828 


0. 


,049 


11 


.80 


2, 


.litO, 


.1 


1.5±0.5 


NGC3221 


0, 


.013 


11 


.00 


0, 


.5±0, 


.2 


<5.9 


NGC3256 


0. 


009 


11 


.56 


66. 


,1±3, 


.1 


<109.3 


ESO264-G036 


0, 


023 


11 


.35 


10. 


.6±1, 


.0 


6.9±3.4 


IRASF10565+2448 


0. 


.043 


12 


.02 


16. 


.6±0, 


.6 


<24.7 


ESO264-G057 


0. 


.017 


11 


.08 


7. 


.3±0, 


.5 


<13.0 


MCG+07-23-019 


0. 


.035 


11 


.61 


0. 


.7±0, 


.1 


<1.5 


CGCGOll-076 


0. 


.025 


11 


.37 


29. 


.9±0. 


.8 


20.0±8.8 


IC2810 


0. 


.034 


11 


.59 


5. 


.9±0, 


.3 


<9.1 


ESO319-G022 


0. 


.016 


11 


.04 


9. 


.4±0, 


.3 


6.4±2.7 


NGC3690/IC694 


0. 


Oil 


11 


.88 


409. 


.0±5. 


.6 


359.1±44.5 


ESO320-G030 


0. 


Oil 


11 


.10 


24. 


.1±1, 


.2 


<35.3 


ESO440-IG058 


0. 


023 


11 


.36 


10. 


,3±0, 


.8 


<21.1 


IRASF12112+0305 


0. 


.073 


12 


.28 


4, 


.6±0, 


.4 


<6.9 


ESO267-G030 


0. 


.018 


11 


.19 


14. 


,1±0, 


.7 


<20.9 


NGC4194 


0. 


.009 


11 


.06 


50. 


.2±3, 


.1 


<77.6 


IRAS12116-5615 


0. 


027 


11 


.59 


42. 


.3±0. 


.8 


29.3±11.6 


IRASF12224-0624 


0. 


026 


11 


.27 


6, 


.2±0, 


.3 


6.1±0.3 


NGC4418 


0. 


.007 


11 


.08 


202.3±12, 


.1 


198.1±12.7 


UGC08058 


0. 


042 


12 


.51 


670. 


,0±2, 


.8 


665.9±4.6 


NGC4922 


0, 


,024 


11 


.32 


43. 


,6±0, 


.5 


39.6±3.6 


CGCG043-099 


0, 


.037 


11 


.62 


9, 


.3±0, 


.4 


<15.1 


MCG-02-33-098/9 


0. 


.016 


11 


.11 


12. 


.0±0. 


.9 


<23.1 


ESO507-G070 


0. 


.022 


11 


.49 


13. 


.6±0. 


.7 


<20.3 


IRAS13052-5711 


0. 


.021 


11 


.34 


7. 


.1±0, 


.4 


<12.5 


NGC5010 


0. 


.021 


11 


.50 


29. 


.2±4. 


.6 


22.9±7.2 


IRAS13120-5453 


0. 


.031 


12 


.26 


45. 


.3±1. 


.4 


<60.2 


IC0860 


0. 


.013 


11 


.17 


6, 


.4±0. 


.2 


4.4±1.8 


VV250a 


0, 


,031 


11 


.74 


14. 


,1±1, 


.1 


<24.9 


UGC08387 


0. 


023 


11 


.67 


19. 


,5±0, 


.9 


<32.7 


NGC5104 


0. 


,019 


11 


.20 


16. 


.1±0, 


.8 


<23.4 


MCG-03-34-064 


0, 


,017 


11 


.24 


211. 


Oil, 


.3 


210.3±1.4 


NGC5135 


0. 


014 


11 


.17 


33. 


,8±1, 


.8 


<45.1 


ESO173-G015 


0, 


,010 


11 


.34 


62. 


,4±2. 


.6 


<85.0 


IC4280 


0. 


016 


11 


.08 


6. 


.5±0, 


.5 


<10.8 


NGC5256 


0. 


028 


11 


.49 


14. 


,0±0. 


.5 


<22.1 


NGC5257/8 


0. 


,023 


11 


.55 


5, 


.3±0, 


.3 


<8.2 


UGC08696 


0, 


,038 


12 


.14 


52. 


,1±1. 


.8 


46.6±5.2 


UGC08739 


0, 


,017 


11 


.08 


9, 


.8±0, 


.5 


<13.5 


ESO221-IG010 


0. 


,010 


11 


.17 


7, 


.7±0, 


.9 


<19.5 


NGC5331 


0. 


,033 


11 


.59 


6. 


.4±0. 


.4 


<11.7 


NGC5394/5 


0. 


.012 


11 


.00 


14. 


.3±1. 


.2 


<23.6 


CGCG247-020 


0. 


.026 


11 


.32 


10. 


.3±0. 


.5 


<16.4 


NGC5653 


0. 


.012 


11 


.06 


7. 


.4±0. 


.5 


<12.2 


IRASF14348-1447 


0. 


.082 


12, 


.30 


7. 


.4±0. 


.4 


5.2±2.0 


IRASF14378-3651 


0. 


.068 


12, 


.15 


6. 


.1±0. 


.4 


<8.0 


NGC5734 


0. 


.014 


11 


.06 


12. 


.8±0, 


.6 


<17.2 


VV340a 


0, 


,033 


11 


.67 


13. 


.4±0, 


.6 


<20.4 


VV705 


0, 


,041 


11 


.89 


1, 


.9±0, 


.3 


<3.0 


ESO099-G004 


0, 


,029 


11 


.67 


15. 


,3±0, 


.5 


<22.9 


IRASF15250+3608 


0, 


,055 


12 


.02 


35. 


,6±2, 


.5 


34.4±2.7 


NGC5936 


0, 


,013 


11 


.07 


11. 


,8±0, 


.9 


<21.6 


UGC09913 


0, 


,018 


12 


.21 


44. 


,7±1, 


.9 


34.9±8.9 


NGC5990 


0, 


,013 


11 


.06 


97. 


,2±0, 


.5 


89.9±6.5 


NGO6052 


0, 


,016 


11 


.02 


3, 


,6±0 


.9 


<10.5 


NGC6090 


0, 


,030 


11 


.51 


1, 


,5±0 


.4 


<8.1 


IRASF16164-0746 


0, 


,027 


11 


.55 


17. 


,3±2, 


.2 


<24.9 


CGCG052-037 


0, 


,024 


11 


.38 


10. 


,5±0, 


.7 


<22.8 


NGC6156 


0. 


Oil 


11 


.07 


16. 


.0±0. 


.3 


12.4±3.2 


ESO069-IG006 


0. 


.046 


11 


.92 


12. 


.4±1. 


.4 


<28.8 
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TABLE 1 

continued - Full sample of infrared galaxies: results of the 

FIT and corrected /gAGN FLUXES 







log(LTIR) 


/gagn 


/eAGN 








mJy 


mJy 


IRASF16399-0937 


0.027 


11.56 


13.5±0.4 


9.7±3.4 


ESO453-G005 


0.021 


11.29 


3.6±0.2 


<5.1 


NGC6240 


0.024 


11.85 


73.5±3.5 


57.4±14.7 


IRASF16516-0948 


0.023 


11.24 


1.7±0.3 


<8.5 


NGC6286 


0.019 


11.32 


14.3±0.6 


<22.1 


IRASF17132+5313 


0.051 


11.89 


3.1±0.2 


<4.8 


IRASF17138-1017 


0.017 


11.42 


5.9±0.9 


<10.5 


IRASF17207-0014 


0.043 


12.39 


19.6±1.1 


<27.4 


ESO138-G027 


0.021 


11.34 


12.2±0.4 


<17.4 


UGC11041 


0.016 


11.04 


8.2±0.6 


<14.3 


CGCG141-034 


0.020 


11.13 


12.5±0.5 


<17.9 


IRAS17578-0400 


0.013 


11.35 


8.9±0.7 


<19.1 


IRAS18090+0130 


0.029 


11.58 


8.5±0.6 


<17.6 


NGC6621 


0.021 


11.23 


8.9±0.6 


<14.6 


CGCG142-034 


0.019 


11.11 


9.6±0.4 


<13.4 


IRASF18293-3413 


0.018 


11.81 


41.1±2.8 


<86.5 


NGC6670A/B 


0.029 


11.60 


9.6±0.6 


<21.3 


IC4734 


0.016 


11.30 


20.2±0.8 


<28.5 


NGC6701 


0.013 


11.05 


16.9±1.0 


<26.0 


ESO593-IG008 


0.049 


11.87 


7.1±0.5 


<12.3 


NGC6786/UGC11415 


0.025 


11.43 


22.1±0.7 


<29.1 


IRASF19297-0406 


0.086 


12.37 


7.7±0.4 


<9.9 


IRAS19542+1110 


0.065 


12.04 


11.1±0.4 


9.3±1.6 


ESO339-G011 


0.019 


11.12 


28.5±0.6 


23.5±4.5 


NGC6907 


0.011 


11.03 


7.5±0.6 


<13.0 


MCG+04-48-002 


0.014 


11.06 


22.6±1.0 


<34.4 


NGC6926 


0.020 


11.26 


5.3±0.1 


4.9±0.4 


IRAS20351+2521 


0.034 


11.54 


5.0±0.5 


<10.7 


CGCG448-020 


0.036 


11.87 


20.2±0.5 


14.3±5.3 


ESO286-IG019 


0.043 


12.00 


39.5±1.0 


37.3±2.2 


ESO286-G035 


0.017 


11.13 


8.1±0.8 


<23.6 


IRAS21101+5810 


0.039 


11.75 


5.6±0.3 


<9.2 


ESO343-IG013 


0.019 


11.07 


12.9±0.5 


<19.7 


NGC7130 


0.016 


11.35 


33.4±0.5 


27.0±5.7 


ESO467-G027 


0.018 


11.02 


3.3±0.3 


<7.0 


ESO602-G025 


0.025 


11.27 


29.9±0.7 


<40.3 


UGC12150 


0.021 


11.29 


17.5±0.8 


<25.7 


ESO239-IG002 


0.043 


11.78 


10.2±0.3 


<13.2 


IRASF22491-1808 


0.077 


12.11 


3.3±0.2 


<4.5 


NGC7469 


0.016 


11.59 


171.8±1.9 


146.7±22.3 


CGCG453-062 


0.025 


11.31 


6.0±0.3 


<9.6 


ESO148-IG002 


0.045 


12.00 


24.7±0.4 


19.1±4.9 


NGC7552 


0.005 


11.03 


114.9±3.5 


<156.9 


IC5298 


0.027 


11.54 


44.1±1.7 


40.2±3.9 


NGC7591 


0.017 


11.05 


12.2±0.5 


<16.5 


NGC7592 


0.024 


11.33 


28.4±0.6 


<37.8 


ESO077-IG014 


0.042 


11.70 


11.2±0.6 


<16.5 


NGC7674 


0.029 


11.50 


138.5±1.2 


137.0±1.8 


NGC7679 


0.017 


11.05 


13.6±1.2 


<24.8 


IRASF23365+3604 


0.064 


12.13 


5.1±0.3 


3.5±1.5 


MCG-01-60-022 


0.023 


11.21 


5.3±0.5 


<12.8 


IRAS23436+5257 


0.034 


11.51 


11.4±0.3 


8.0±3.0 


NGC7752/3 


0.017 


11.01 


5.3±0.4 


<9.5 


NGC7771 


0.014 


11.34 


17.8±1.0 


<26.8 


MRK0331 


0.018 


11.41 


22.3±1.5 


<46.9 



Sources from lTran et al.1 (120011 ) 



IRASF00183-7111 


0, 


,327 


12, 


.77 


45, 


,1±0. 


,6 


44.9±0.6 


IRAS00188-0856 


0, 


129 


12, 


.31 


11, 


,8±0. 


,3 


11.1±0.7 


IRAS00275-2859 


0, 


,279 


12, 


.46 


35, 


,8±0. 


1 


35.4±0.3 


IRAS00406-3127 


0, 


,342 


12, 


.64 


11, 


,6±0. 


3 


11.5±0.3 


IRAS02113-2937 


0, 


194 


12, 


.29 


1, 


.7±0, 


1 


<2.3 


IRASF02115+0226 


0, 


,400 


12, 


.48 


1, 


,2±0. 


1 


1.0±0.2 


IRASF02455-2220 


0, 


,296 


12, 


.57 


1, 


,8±0, 


1 


1.6±0.2 


IRAS03000-2719 


0, 


,221 


12, 


.41 


2, 


.7±0, 


1 


2.2±0.5 


IRAS03538-6432 


0, 


,310 


12, 


.65 


8, 


,9±0, 


,2 


8.4±0.5 


IRAS03521+0028 


0, 


152 


12, 


.46 


2, 


,1±0. 


1 


1.5±0.6 


IRAS04384-4848 


0, 


,213 


12, 


.32 


3, 


,8±0. 


,3 


3.3±0.5 


IRAS17463+5806 


0, 


,309 


12, 


.48 


2, 


.2±0, 


1 


1.9±0.3 


IRAS18030+0705 


0, 


146 


12, 


.18 


0, 


,5±0. 


1 


<2.7 


IRASF23529-2119 


0, 


,430 


12, 


.55 


8, 


.9±0, 


1 


8.8±0.1 
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TABLE 2 

Distributions parameters of measurements and simulations 



log(LTIR) 
Lq 


Measurements 


Detect. limit 
mJy 




Simulations 


m 


cr 




0"dct 


d.r.|7oJ 




""det 


d.r.|%J 


10.9-M1.3 


-1.9 


0.6 


23.3 


17.9 


-2.0 


0.6 


23.4 


-3.1 


1.0 


11.3-M1.7 


-1.6 


0.5 


26.8 


17.5 


-1.9 


0.6 


26.9 


-3.0 


0.9 


11.7-M2.1 


-1.8 


0.5 


48.0 


10.5 


-2.0 


0.6 


48.2 


-2.6 


0.8 


12.1-M2.5 


-1.7 


0.6 


65.2 


8.3 


-1.9 


0.7 


65.2 


-2.3 


0.9 


12.5-M2.9 


-1.0 


0.5 


100.0 


0.9 


-1.0 


0.5 


100.0 


-1.0 


0.5 



TABLE 3 

High redshift sample of infrared galaxies: results of the fit and 
corrected /bagn fluxes 



Name 




2 


log(LTIR) 


/SAGN 


/6.2pcakSB 


/gagn 










mJy 


mJy 


mJy 






Sources fromlBrand et al.l (l2008l) 




70Bootesl 




0.50 


12.15 


0.21±0.04 


0.30±0.08 


0.18±0.05 


70Bootes2 




0.37 


11.79 


0.12±0.04 


0.50±0.07 


<0.12 


70Bootes3 




0.99 


13.05 


0.15±0.04 


1.10±0.08 


<0.25 


70Bootes4 




0.98 


12.95 


0.14±0.04 


1.15±0.07 


<0.14 


70Bootes5 




1.21 


13.24 


1.98±0.08 


0.28±0.14 


1.96±0.08 


70Bootes6 




0.94 


13.01 


1.81±0.07 


0.25±0.12 


1.79±0.07 


70Bootes7 




0.66 


12.61 


0.32±0.04 


1.30±0.07 


<0.43 


70Bootes8 




0.81 


12.88 


0.74±0.06 


0.30±0.10 


0.71±0.06 


70Bootes9 




0.67 


12.81 


0.23±0.05 


2.38±0.09 


<0.42 


70Bootesl0 




0.51 


12.41 


0.72±0.06 


0.29±0.10 


0.69±0.07 


70Bootesll 




0.48 


12.49 


2.01±0.07 


0.31±0.13 


1.98±0.08 






Sources from iPoDe et all (1200811 




C3 




1.88 


12.78 


0.15±0.03 


0.10±0.05 


0.14±0.03 


GN39 




1.98 


12.70 


0.08±0.04 


0.54±0.07 


<0.08 


GN07 




1.99 


12.84 


0.09±0.03 


0.30±0.06 


<0.09 


GN06 




2.00 


12.81 


0.06±0.03 


0.51±0.05 


<0.06 


CI 




2.01 


12.98 


0.48±0.07 


0.30±0.12 


0.45±0.07 


GN26 




1.23 


12.56 


0.05±0.02 


0.81±0.03 


<0.14 


GN17 




1.73 


12.30 


0.17±0.03 


0.39±0.06 


0.13±0.05 


GN19 




2.48 


13.08 


0.13±0.04 


0.26±0.06 


<0.13 


GN04 




2.55 


12.82 


0.24±0.04 


0.15±0.08 


0.23±0.05 




Sources from iMenendez-Delmestre et al.l (120071) 




SMMJ030228- 


f000654 1.41 


13.44 


0.07±0.03 


0.27±0.05 


<0.12 


SMMJ163639- 


f 405636 1.50 


12.81 


0.09±0.02 


0.17±0.04 


0.07±0.03 


SMMJ163650- 


f405735 2.38 


13.52 


0.41±0.05 


0.33±0.07 


0.38±0.06 






Sources from IValiante et al.l (120071) 




SMMJ00266-I-1708 


2.73 


12.70 


0.15±0.05 


0.30±0.09 


<0.15 


SMMJ02399-0136 


2.81 


13.40 


0.42±0.04 


0.21±0.07 


0.40±0.04 


SMMJ09429-I-4659 


2.38 


13.80 


0.08±0.04 


0.43±0.07 


<0.08 


SMMJ09431-I-4700 


3.36 


12.95 


0.58±0.04 


0.29±0.07 


0.55±0.04 


SMMJ10519-f5723 


2.67 


13.07 


0.11±0.03 


0.14±0.05 


0.10±0.03 


SMMJ1052H-5719 


2.69 


13.25 


0.13±0.03 


0.13±0.05 


0.11±0.03 


MMJ154127-I-6616 


2.79 


12.95 


0.14±0.01 


0.19±0.02 


0.12±0.02 


SMMJ16369-I-4057 


1.21 


12.23 


0.15±0.02 


0.23±0.03 


0.12±0.03 


SMMJ16371-I-4053 


2.38 


13.04 


0.14±0.03 


0.18±0.05 


0.12±0.03 



TABLE 4 
Evolution models 



Model 


Description 


rii 


n2 


21 


mi 


m2 


22 


Ml 


SB 


3.6 


0.0 


2.0 


1.0 


0.0 


2.0 


M2 


SB-I-AGN 


3.6 


0.0 


2.0 


1.0 


0.0 


2.0 


M3 


SB-I-AGN 


3.4 


-1.0 


2.3 


2.0 


-1.5 


1.0 


M4 


SB-I-AGN-I- L-T evolution 


3.4 


-1.0 


2.3 


2.0 


-1.5 


1.0 


M5 


SB-I-AGN -1- L-T evolution+ "cold spread" 


3.4 


-1.0 


2.3 


2.0 


-1.5 


1.0 


M6 


SB-|-"cut" AGN-I-L - T evolution+ "cold spread" 


3.4 


-1.0 


2.3 


2.0 


-1.5 


1.0 
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Fig. 1. — Total infrared luminosit y distribution of the galaxies of our sample, including galaxies from RBGS catalog l|Sanders et al.ll2003l ) 
and sources from lTran et al.l 1 I2OOIII . 
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Rest wavelength [^m] Rest wavelength [/im] 




0.05 - - 0.2 - 



0.00 L ^ „ I . ^ . „ ^ „ . . i no t , I . , . . ' 

5.8 6.0 6.2 6.4 6.6 6.8 5.8 6.0 6.2 6.4 6.6 6.8 

Rest wovelength [jum] Rest wovelength [^m] 

Fig. 2. — Examples for the decomposition used to isolate the AGN continuum. Top panels: Spitzer IRS low resolution spectra of four 
representative galaxies. Bottom panels: Cutout of the region around the 6.2 jim PAH feature. Top continuous line = observed spectrum. 
Top dotted line = fit by the sum of the M82 spectrum and a linear AGN continuum. Bottom dotted line = fitted AGN continuum. Bottom 
continuous line = difference of observed spectrum and fitted PAH component. The thick vertical line indicates the 6 /um AGN continuum 
flux density /eAGN- 
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log(i'Lw'L,.) \oq{vLjU) 



12.5<log(L™ 


)<1 


2.9 - 


m«=-1-0 
<J,=OJ 






<l.r..lOO.OX 















-4 -3 -2 -1 

log(vL»/L™) 



Fig. 3. Distribution of i^ieAGN/i'TlRi for different luminosity bins. Solid histograms show the measurements, while the cross-hatched 
histograms indicate objects with upper limits. As indicated in each diagram, the mean values of the detections do not vary significantly 
between the different bins. What noticeably varies is the detection rate, spanning values from 23.3% to 100%. 
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log(<.U/L,.) log(i.U/L,.) 



2500 



2000 
1500 
1000 
500 





m = -1.0 
a=0.5 



o.=0.5 
().r,= 100.0X 



I 



-4 -2 

log(i.L^L™) 



Fig. 4. — Simulations of i/LgAGN/^'TlRi assumed to be normally distributed, reproducing the same detection rates (d.r^, means (rrtdct) 
and sigma (iTdct) of the detections as the data. The detection limits assumed for /eAGN in each bin is reported in Tab. [T] Symbols have 
the same meaning as in Fig. [3] Each simulation is made for a population of 2 X 10^ objects. Each diagram reproduces the values observed 
in the corresponding panel of Fig. [3] and shows also mean (m) and sigma (cr) of the adopted distribution. 
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1 1.7<log(UR)<12.5 


ni.= -1.3 

0„=0.4 

d.r.=83.3X 






















m, 








4 -3 


-2 




-1 


C 



log(i/L4./L„.) 




-4 -3 -2 -1 

109(1/ L»./L™) 



13.3<log(UR)<14.1 


m,= -I,l 
a.=0.2 


m 








d.r.=50.05! 




















4 -3 


-2 




-1 


c 



log(i/L»yL,,) 



Fig. 5. — Distribution of i^igAGN/^TlR for different luminosity bins. The sample lies at redshift 0.37 < z < 3.35. Symbols have the 
same meaning as m Fig. [3] Each dia gram also shows the mean, the dispersion and the detection rate of the distribution. 
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Fig. 6. — Simulations of J^isAGN/^TlR normal distributions reproducing the same detection rates (d.r.), means (rridct) a^nd sigma (fJ^ct) 
of the detections obtained from the data of the high redshift sample, assuming a detection limit of O.lSmJy. Each diagram reproduces the 
values observed in the corresponding panel of Fig. [5] while the last luminosity bin has not been simulated because the statistic is too low 
for the fit to converge. Symbols have the same meaning as in Fig. [3] Each diagram also shows mean (m) and sigma (cr) of the adopted 
distribution. 
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Fig. 7. — J^LgAGN/i'TlRi vs. Ltir- Dots and small arrows represent measurements and upper limits of all the galaxies in our local 
sample. The large open squares are the means of the detections (m^pt) in bins of luminosity (see Tab. [T] and Figl^l. Due to the large 
number of non-detections, they represent only an upper-limit on the true log i^L^agn/^tir for almost all the bins. The green stars and the 
associated bars show means (m) and sigma (cr) of the gaussian distributions adopted in the simulations in order to reproduce detection rates 
(d.t.), means (mdot) and sigma (crdct) of the detections from the measurements (see Tab.[l]and Fig|4)l. The best fit, i'Lsagn/^tir i'xiR' 
gives a = 1.4 ± 0.6 and is shown as a dotted green line: the relation is assumed to be flat at high Ltir- The red stars and their bars show 
means and sigma of the gaussian distributions adopted to reproduce measurements of the high redshifts sample (see Fig. [5] and Fig. |6j. 
The relation found for the local galaxies no longer holds. 
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Fig. 8. — Relation between the color ij(60, 100) and -Ltir (solid and dotted lines give the mean and 1— cr envelope of Eq. [2] respectively). 
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Fig. 9. — Some ex amples of adopted SEDs (solid lines) ob tained adding a star-fo rming galaxy spectrum (dashed lines) IIDale et al.ll200ll : 
IDale fc Helou ||2002I) and an AGN template (dotted lines) (lEfstathiou et aT1ll995l) . The flux density is expressed in Jansky and is scaled 
for a distance of 100 Mpc. In each diagram we indicate Ltir E'nd the ratio between the luminosity of the AGN and the luminosity of 
the starburst spectra. For each luminosity class, two examples are shown: an AGN-dominated SED (iAGN/^SB ^ li left column) and a 
starburst-dominated SED (I/AGN/i'SB < Ij right column). 
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Fig. 10. — The 'double power-law" parametrization of the TIR luminosity function. 
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Fig. 11. — Number counts at 24, 70, 160, 850 /im and redshift distributions at 24 and 850 fi m, considering evoluti ons with ni = 3.6, 
712 = 0, mi = 1, 971 2 = 0, 21 = 22 = 2, without AGN contribution (Ml, solid lin es). Data are from lPapovich et al.l ||2004|') (filled circles) and 
IShupe et all ll200l) ( open squares) fo r the 24 /im number co unts. [Fraver et al.| ir2 006a b) (70 fim and 160 fim counts), ICoppirTe t al. ( 200(|) 
(850 /.tm counts'). IWuvts et al.l l|2008 ) [24 fira redshifts) and IChapmanet'arr i 20051 ) (850 ^im redshifts). EucUdean normalized differential 
counts are shown at 24, 70, 160 fira, while integrated counts are shown at 850 /im following the style of the papers in which they were 
originally presented. Redshift data (open histograms) have been scaled for comparison with the simulations {filled histograms). 
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Fig. 12. — Same as Fig. 1111 but now with AGN contribution (M2). AGN-dominated galaxies (Lagn/^'SB > 1) are represented with 
dashed line, while starburst dominated galaxies (iAGN/^SB < 1) are in dotted line. 
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Fig. 13. — Number counts at 24, 70, 160, 850 /im and redshift distributions at 24 and 850 fim, consid erin g evolutions with ni = 3 .4, 
n2 = —1, 21 = 2.3, mi = 2, m2 = —1.5 and Z2 = 1, with AGN contribution (M3). Symbols arc like in Fig. [12] Data arc like in Fig. [iT] 
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Fig. 14.— 



Same as Fig. 1131 assuming an evolution in the L — T relation (M4). 
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Fig. 15.— 



Same as Fig. 1131 assuming an evolution and a spread towards lower temperatures in the L — T relation (M5). 
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Fig. 16. — Cosmic background from mid-IR to millimeter wavelengths for model M5. The CIB derived from t he models at 24, 70 , 
100, 160 and 850 ^m is shown by big open squares. The analytic form of the CIB at the FIRAS wavelengths is from'Fixsen et aTl I I1998I ) 
(solid an d dashed line s). Measure ments are at DIRBE wavebands (140 and 240 /^m, plus si gns, Lagache et al. 2000), at 100 /^m [triangles, 
ILagache et al.ii2000i ; [Renault et al.ll2001ll and at 60 fim (frmnq/eg. IMiville-Deschenes et al.ll2002ll . See iLagache et al.i (I2003ii'l for references 
about the limits at shorter wavelengths. 
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Fi g. 17. — Luminosi ty function of model M5 calculat ed at differe nt redshifts, compare d with av ailable measu r ement s. Data are 
from IHuang et al.l 1120071') (squares), ILe Floc'h et al.l II2005I ') (triangles), IMagnelli et al.l II2009I ') (circles), ICaputi et"an II2007I') (diamonds), 
IChapman et al.l II2005I) (stars) . 
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Fig . 18. — Histogram showing the SFR evolution of the model M5 compared with the observational results shown bv IChapman et al.l 
II2005I ) (see their Fig. 12). Yellow triangle: radio surveys; red circles: mid-IR surveys ; green stars UV and optical surveys, corrected for 
dust extinction; violet squares: SMGs corrected for completeness. See IChapman et al.l (120051 '! for a full list of references. 
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Fig. 19. — Fraction of all sources whose mid-IR emission (a) and total infrared emission (6) are dominated by AGN as a function of /24 
for model M5 (sol i d hist ogram) and M6 [dashed histogram). Redshift range of the sources is < 2: < 8. Open squares represent the results 
from lBrand et al.l igOOi). 
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Fig. 20. — Predictions using the best-fit model M5 for the upcoming Herschel surveys. Top: integral number counts at 100 /im {solid 
line), 250 /im {dotted line), 350 /im {dashed line), 500 /tm {dot-dashed line). Middle: redshift distributions of PEP surveys in the deep 
(10' X 15', /loo > 1-5 mjy, left) and large (85' X 85', /leo > 8mjy, right) fields. Bottom: expected distributions for the survey in the CDFS 
field for model M4 {dotted line) and M5 {solid line) of redshifts {left) and temperatures {right). 



